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1.  Scope  of  the  Program 

The  program  is  a  general  purpose  finite-difference  FORTRAN-77 
code  designed  to  solve  steady-state  three-dimensional  heat  transfer  and 
fluid  flow  problems  in  rectangular  coordinates.  The  code  is  designed  to 
effectively  handle  conjugate  domain  problems,  arising  when  both  solid  and 
fluid  are  present  in  the  computational  region.  Pure  heat  conduction 
problems  or  only  fluid  flow  problems  can  also  be  solved. 

A  control-volume  formulation  is  used  to  discretize  the  governing 
equations  expressed  in  a  canonical  form.  The  control  volumes  for  the 
velocity  components  are  staggered  with  respect  to  those  for  the  temperature 
and  pressure.  Harmonic  mean  formulation  for  the  interface  diffusivities  is 
used  to  effectively  handle  sharp  discontinuities  in  property  values  in  the 
computation  domain,  without  having  to  resort  to  very  fine  mesh  in  the 
region  of  discontinuities.  The  velocity-pressure  coupling  is  handled  by  the 
SIMPLER  algorithm  of  Patankar  [1],  The  solution  to  the  steady  state 
problem  is  achieved  by  providing  an  initial  guess  of  all  the  dependent 
variables  and  proceeding  through  an  iterative  scheme  using  the  Thomas 
algorithm  (also  known  as  the  Tri-Diagonal  Matrix  Algorithm).  Iterations 
are  to  be  stopped  when  values  of  variables  in  successive  iterations  do  not 
change  more  than  a  specified  convergence  criterion. 

2.  Program  Structure 

The  code  is  highly  modular  consisting  of  a  short  main  program  and 
various  subroutines.  The  problem  independent  portion  consists  of  the  main 
program  and  subroutines  PROFIL,  TDMA,  FORM  and  OPTION.  The 
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problem  dependent  part,  which  the  user  is  primarily  concerned  with, 

consists  of  a  BLOCK  DATA  subprogram  and  a  subroutine  USE. 

A  brief  description  of  the  various  subprograms/subroutines  follows: 

PROFEL  The  profile  assumptions  in  between  node-points  for  the 

dependent  variables  are  incorporated  here.  Presently,  it 
contains  the  power  law  profiles  [1].  Other  profile  schemes 
such  as  the  hybrid  scheme  or  the  central  differencing 
scheme  may  be  incorporated  by  the  user. 

TDMA  The  line-by-line  TDMA  (Tri-Diagonal  Matrix  Algorithm)  is 

used  for  solving  the  equations.  The  subroutine  solves  for  a 
canonical  equation  described  later  (see  Eq.  1). 

FORM  Evaluates  the  geometric  parameters  such  as  lengths,  areas 

etc.  in  ENTRY  GEOMET  and  the  coefficients  for  the  various 
equations  in  ENTRY  COEFF  which  are  then  used  to  solve 
for  various  equations  in  accordance  with  the  SIMPLER 
algorithm  [1]. 

OPTION  The  user  may  access  two  optional  utilities  (entries), 

UMESH  and  PRINT.  UMESH  is  to  be  accessed  if  a  uniform 
grid  is  required,  while  PRINT  is  called  if  dependent 
variables  are  to  be  printed  out  over  the  entire  domain. 
Calling  PRINT  causes  a  plane  by  plane  field  printout,  for 
different  planes  perpendicular  to  the  z  axis. 

BLOCK  DATA  Contains  block  COMMONS  that  contain  default  values  of 
various  quantities  described  subsequently.  The  default 
values  can  be  changed  by  the  user. 
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USE 


Contains  problem  dependent  information  to  be  specified  by 
the  user.  The  different  entries  in  the  subroutine  USE  are 
MESH,  BEGIN,  VARRHO,  BNDRY,  PRTOUT  and  DIFFUS 
The  usage  of  these  entries  is  described  later. 

3.  Program  Usage 


It  is  required  at  first  to  cast  the  equations  to  be  solved  in  the  following 
canonical  form: 


^(Pu<t>)  +  ^(PV(t>)  +  ^(pw<j>) 


— (r— +  -2-(r— i  +  —  (r— )  +  s 

dxU  Ox  +  3yu  8y'  +  0z(1  3z' 


(i) 


where  p  is  the  fluid  density,  u,  v  and  w  are  the  velocity  components  in  the  x, 
y  and  z  directions  respectively,  <j>  is  the  dependent  variable  to  be  solved  for,  F 
is  the  appropriate  diflusivity  and  S  is  the  source  term. 

For  efficient  execution  of  the  algorithm,  and  for  better  convergence 
prospects,  it  is  necessary  that  S  be  represented  in  the  form 

S=Sc+Sp<j)  (2) 


and  it  is  required  that  Sp  <  0.  For  an  appropriate  linearization  of  the 
source  term  S,  when  it  does  not  have  a  readily  obtainable  form  of  Eq.  (2),  the 
user  is  referred  to  Patankar  [1],  pp.  48-49. 
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It  is  not  required  to  cast  the  continuity  equation  into  the  canonical 
form,  since  this  is  incorporated  in  the  program  internally.  For  example, 
when  solving  a  combined  incompressible  fluid-flow,  heat  transfer  problem, 
the  user  needs  to  cast  the  three  momentum  equations  and  one  energy 
equation  in  the  canonical  form  of  Eq.  (1).  Note  that  the  pressure  gradient 
terms  encountered  in  the  momentum  equations  are  incorporated  internally 
and  the  user  is  cautioned  not  to  incorporate  them  in  the  source  term  S. 

As  mentioned  earlier,  the  user  needs  to  be  concerned  only  with  the 
BLOCK  DATA  and  subroutine  USE.  A  detailed  description  of  these 
segments  now  follows. 


BLOCK  DATA 

The  array  F(I,J,K,NF)  contains  different  variables  or  properties.  The 
default  equivalence  of  the  array  F  is  given  below: 


NF  Quantity  Stored  in  F(I.J.K.NF)  FORTRAN  Array  in  USE 


1  velocity  in  x  direction,  u 

2  velocity  in  y  direction,  v 

3  velocity  in  z  direction,  w 

4  pressure  correction 

5  temperature 

6  pressure 


U(I,J,K) 

V(I,J,K) 

W(I,J,K) 

PC(I,J,K) 

T(I,J,K) 

P(I,J,K) 
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The  BLOCK  DATA  contains  several  data  that  may  be  reinitialized  by 

the  user.  Each  of  these  is  explained  below: 

IPREF,  JPREF,  KPREF  - 

Values  of  the  indices  (I,J,K)  where  the  pressure  is 
set  to  zero.  Note  that  the  program  evaluates  only 
pressure  differences  and  not  the  actual  pressures. 

NTIMES  -  The  value  of  NTIMES(NF)  determines  the  number 

of  sweeps  of  the  TDMA  solution  over  the  entire 
domain  during  one  iteration.  More  than  one 
sweep  per  iteration  may  be  required  for  highly 
nonlinear  equations  to  assure  convergence. 
Clearly,  a  larger  value  of  NTIMES(NF)  will  result 
in  a  longer  computational  time. 

LAST  -  Defines  the  total  number  of  iterations  to  be 

performed. 

LPRINT  -  The  default  value  of  this  logical  FORTRAN 

variable  is  FALSE.  Setting  LPRINT(NF)  to  be  true 
will  cause  a  field  printout  of  the  particular 
variable  associated  with  NF  when  a  call  to  the 
entry  PRINT  is  invoked  through  the  statement 
CALL  PRINT. 

LSOLVE  -  The  default  value  of  this  logical  FORTRAN 

variable  is  FALSE.  Setting  LSOLVE(NF)  to  be  true 
causes  the  solution  of  the  particular  variable 
associated  with  NF.  In  a  pure  conduction 
problem,  for  example,  only  LSOLVE(5)  must  be  set 
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RELAX - 


RHOCON - 
TITLE(NF)  - 


true.  In  a  pure  fluid  flow  problem  (no  heat 
transfer)  LSOLVE(NF),  NF=1,2,3,4,  and  6  must  be 
declared  true.  In  a  heat  transfer  and  fluid  flow 
problem,  LSOLVE(NF),  NF=1  to  6  should  be  set 
true. 

The  value  for  RELAX(NF)  is  the  underrelaxation 
factor  for  the  solution  of  the  variable  associated 
with  NF.  RELAX(NF)>1  causes  overrelaxation, 
while  RELAX(NF)<1  causes  underrelaxation. 
Underrelaxation  may  be  necessary  in  inany 
problems  when  convergence  is  difficult  to  obtain, 
and  will  increase  the  computational  time. 

Value  of  p  in  Eq.(l)  is  specified  here. 

This  describes  the  title  to  be  given  to  each  field 
printout. 


Subroutine  USE  contains  various  segments  (entries)  called  by  the 
main  program.  The  purpose  of  these  and  their  usage  is  described  here. 

See  Appendix  A  for  meaning  of  the  various  symbols. 

MESH:  The  user  supplies  a  mesh  in  this  entry.  The  computational 

domain  must  be  divided  into  rectangular  control  volumes.  The 
internal  grid  points  are  placed  at  the  geometric  centers  of  the 
control  volumes.  The  boundary  grid  points  are  placed  on  the 
geometric  centers  on  the  control  volumes  faces  that  form  part 
of  the  boundary  (See  Figs.  A-l  and  A-2  in  Appendix  A). 
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Two  options  are  available- 

1.  Uniform  grid  -  The  user  supplies  values  of  LI,  Ml  and  N 1 ; 
XL,  YL  and  ZL  (see  Appendix  A);  and  then  invokes  CALL 
UMESH  to  form  a  uniform  mesh,  i.e.  all  the  control  volumes  in 
a  particular  direction  have  the  same  length: 

ENTRY  MESH 

Ll=10 

Ml=10 

Nl=20 

XL=2. 

YL=3. 

ZL=1. 

CALL  UMESH 
RETURN 

This  implies  a  computational  domain  that  is  2,  3  and  1  units 
long  in  the  x,  y  and  z  directions  respectively.  All  the  control 
volumes  thus  have  dimensions  of  2/10,  3/10  and  1/20  units  in 
the  x,  y  and  z  directions  respectively. 

2.  Non-uniform  grid  -  The  user  supplies  values  of  (XU(I), 

1=1, LI),  (YV(J),  J=1,M1),  ZW(K),  K=1,N1);  XL,  YL,  ZL;  LI,  Ml 
and  Nl.  The  code  internally  determines  the  locations  of  the 
boundary  and  internal  grid  points  (stored  in  X(I),Y(J),Z(K)) 
which  are  placed  at  the  geometric  centers  of  the  control 
volumes.  For  example,  X(I)=0.5*[XU(I)+XU(I+1)]  etc. 
Boundary  points  are  grid  points  that  are  located  on  planes  x=0 
or  x=XL  or  y=0  or  y=YL  or  z=0  or  z=ZL. 

BEGIN:  The  user  specifies/calculates  quantities  that  are  required  prior 

to  the  start  of  iterations  such  as  parameter  values,  initial 
guess  to  the  solution  and  boundary  conditions. 
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Note: 

VARRIIO: 

BNDRY: 

Note: 


Entries  MESH  and  BEGIN  are  called  by  the  main  program 
only  once  before  the  iterations  begin. 


Functional  dependency  of  the  fluid  density  can  be  accounted  for 
in  this  entry  for  variable  density  situations. 


The  boundary  values  of  the  variables  are  updated  every 
iteration  in  accordance  with  the  boundary  conditions. 

Boundary  conditions  that  involve  specification  of  the  dependent 
variable  value  can  be  incorporated  in  BEGIN  and  it  is  not 
necessary  to  include  these  in  BNDRY.  This  is  so  because  the 
program  solves  for  the  dependent  variables  only  at  the  internal 
grid  points  and  not  at  the  boundary  points. 

Examples  of  different  boundary  conditions: 


ENTRY  BNDRY 

DO  1 1=1, LI 
DO  1  J=1,M1 

C  Insulated  boundary 

C  At  z=ZL,  dT/dz=0 

T(I,J,N  1)=T(I,  J,N  1- 1 ) 

C  Convective  condition 

C  At  z=0,  h(Tf-T)=-kdT/dz 

C  H=h  (convection  coefficient),  TF=Tf  (fluid  temperature) 

C  TK=k  (thermal  conductivity) 

DELZ=Z(2)-Z(1) 

T(I,J,1)=(H*DELZ*TF+TK*T(I,J,2))/(TK+H*DELZ) 

C  Specified  heat  flux  condition 

C  at  z=ZL,  -kdT/dz=-q”,  q">0 

C  QPRIME=q"  (heat  flux) 

DELZ=Z(N  1)-Z(N  1-1) 

T(I,J,N1)=T(I,J,N1-1)+QPRIME*DELZ/TK 
1  CONTINUE 


R 


RETURN 


PRTOUT: 


DIFFUS: 


It  may  be  necessary  to  monitor  values  of  different  variables  or 
other  quantities  of  interest  as  each  iteration  proceeds.  This 
can  be  done  in  ENTRY  PRTOUT.  At  the  end  of  the  iterations, 
CALL  PRINT  invokes  a  field  printout  of  the  variables  for  which 
LPRINT  has  been  set  true. 


The  user  specifies  T,  Sc  and  Sp  (Eqs.  1,2)  in  the  variables 

GAM(I,J,K),  CON(I,J,K)  and  AP(I,J,K),  respectively.  The 
diffusivity  T  is  specified  at  all  points,  while  Sc  and  Sp  are 

specified  at  the  internal  points  (i.e.  all  points  excluding  the 
boundary  points).  For  example: 


ENTRY  DIFFUS 

C  Specify  NF  since  DIFFUS  is  called  by  the  main  program 
C  more  than  once  every  iteration. 

DO  2  1=1, LI 
DO  2  J=1,M1 
DO  2  K=1,N1 

C  VISC  is  the  viscosity 

IF((NF.EQ.l).OR.(NF.EQ.2).OR.(NF.EQ.3))GAM(I,J,K)= 

1VISC 

C  TK  is  the  thermal  conductivity,  CP  is  the  specific  heat 
C  capacity. 

IF(NF.EQ.5)GAM(I,J,K)=TK/CP 
2  CONTINUE 

ENDIF 

C  In  the  event  of  uniform  volumetric  heat  generation  in 
C  the  domain  (Q  W/m^). 

DO  3  I=2,L1-1 
DO  3  J=2,M1-1 
DO  3  K=2,N1-1 

IF(NF.EQ.5)CON(I,J,K)=Q/CP 
C  If  AP(I,J,K)  is  not  specified,  it  is  zero  by  FORTRAN 
C  default. 
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CONTINUE 


C  .  SPECIAL  TREATMENT  OF  FLUX  TYPE  BOUNDARY 

C  CONDITIONS  (TO  IMPROVE  CONVERGENCE) 

C  SET  BOUNDARY  VALUE  OF  DIFFUSIVITY  TO  ZERO. 
C  CONVERT  SURFACE  HEAT  FLUX  TO  AN 

C  EQUIVALENT  HEAT  SOURCE  IN  THE  CONTROL 

C  VOLUME  ADJACENT  TO  THE  BOUNDARY. 

DO  4  1=1, LI 
DO  4  J=1,M1 

C  Insulated  boundary 

C  At  z=ZL,  dT/dz=0, 

IF(NF.EQ.5)GAM(I,J,N1)=0.0 
C  Convective  condition 

C  At  z=0,  h(Tf-T)=-kdT/dz 

C  S=Equivalent  heat  source/volume=surface  heat  flux* 

C  surface  area/volume  (source  term  for  next  to  boundary 
C  node). 

C  Surface  convective  flux=h(Tf-T)=(Tf-Ti)/[l/h+(zj-zb)/k] 

C  where  T{  is  the  temperature  at  the  node  adjacent 

C  to  the  boundary,  zj  is  the  z  location  of  the  node,  zb  is  the 
C  z  location  of  the  boundary.  The  factor  in  square  brackets 
C  is  the  thermal  resistance  (Rth)  due  t°  convection  and 

C  conduction. 

C  Surface  area=XCV(I)*YCV(J). 

C  Volume  of  control  volume=XCV(I)*YCV(J)*ZCV(2). 

C  Surface  area/volume=l./ZCV(2) 

C  Tj  =  T(I,J,2),  Zi-zb=Z(2)-Z(l),  T=T(I,J,1). 

C  Source  S  may  be  split  into  Sc  and  Sp 

C  Sc=(Tf/Rth)/ZCV(2),  SD=-l/(Rth*ZCV(2)) 

RTH=(  l./H+(Z(2)-Z(  1))/TK) 

IF(NF.EQ.5)THEN 

GAM(I,J,1)=0.0 

C  Add  Sc  and  Sp  for  the  control  volume  adjacent  to  the 

C  boundary. 

CON(I,J,2)=CON(I,J,2)+TF/(RTH*ZCV(2)) 

AP(I,J,2)=AP(I,J,2)-1./(RTH*ZCV(2)) 

ENDIF 

C  Specified  heat  flux  condition 

C  At  z=ZL,  -kdT/dz=-q",  q">0 

C  QPRIME=q"  (heat  flux) 

IF(NF.EQ.5)THEN 

GAM(I,J,N1)=0.0 

CON(I,  J,N  1-1  )=CON(I,  J,N  1  - 1  )+QPRIME/ZCV(N  1- 1 ) 
ENDIF 

4  CONTINUE 

RETURN 


to 


Conjugate  Domains  Frequently  it  is  required  to  calculate  heat  transfer 
in  domains  that  contain  both  fluid  and  solid.  The  solid  then  can  be  easily 
modeled  by  setting  the  fluid  viscosity  to  a  very  large  value  (eg.  10^0)  in  the 
solid  region  of  the  computational  domain.  The  velocities  in  that  region  then 
will  be  computed  by  the  code  to  be  zero.  Appropriate  solid  properties  such 
as  thermal  conductivity  can  be  prescribed  in  the  solid  region.  Care  must  be 
taken  that  the  control  volume  faces  coincide  with  the  actual  physical 
discontinuities  in  the  domain. 

4.  Program  Flow-chart 

The  flow-chart  for  the  program  execution  is  shown  in  Fig.  1.  The 
usage  and  the  purpose  of  the  various  subroutines  and  subprograms  has 
already  been  explained  in  the  preceding  sections.  A  single  iteration 
consists  of  execution  of  the  SIMPLER  algorithm  to  evaluate  the  velocity  and 
pressure  field  in  the  domain  and  then  the  solution  of  any  other  dependent 
variables  to  be  solved  for,  based  on  the  calculated  velocities. 


li 


Fig.  1  Program  Flow-chart. 
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APPENDIX  A  -  Other  FORTRAN  Variables  in  USE 
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ITER 


Counter  for  iterations.  Program  stops  when 
ITER=LAST. 

LI, Ml, NT  Number  of  grid  points  in  the  X,  Y  and  Z  directions, 

respectively.  Ll-2,  Ml-2  and  Nl-2  are  the  total 
number  of  control  volumes  in  the  X,  Y  and  Z 
directions. 

XL,YL,ZL  Domain  lengths  in  the  X,Y,  and  Z  directions, 

respectively. 

X(I),Y(J),Z(K)  Locations  of  the  node  points  where  temperatures  and 

pressures  are  stored.  The  temperature  T(I,J,K)  is 
stored  at  the  location  (X(I),Y(J),Z(K)).  Note 
X(1)=0.,X(L1)=XL. 

XU(I),YV(J),ZW(K)  Locations  of  control  volume  faces.  XU(I)  is  the 

location  of  the  control  volume  face  perpendicular  to 
the  X  axis  for  the  control  volume  around  the  point 
(X(I),Y(J),Z(K))  such  that  XU(I)<X(I). 

XU(2)=0.,XU(L1)=XL.  XU(1)  is  meaningless. 

U(I,J,K)  is  evaluated  at  the  point  (XU(I),Y(J),Z(K)). 

Similar  description  is  applicable  for  YV(J)  and 
ZW(K). 

FX(I),  FXMfl)  Interpolation  factors  in  the  x  direction  that  enable 

calculation  of  a  quantity  a  control  volume  interface 
perpendicular  to  the  x  direction.  EXAMPLE:  | 

Temperature  at  the  interface  of  two  control  volumes, 
say  at  XU(I),  can  be  calculated  as 

FX(D*T(I,J,K)+FXM(I)*T(I-1,J,K) 
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XDIF(I),YDIF(J), 

ZDIF(K) 

XCV(I),  YCV(J), 
ZCV(K) 


Similar  interpretation  for  FY(J),  FYM(J);  FZ(K), 
FZM(K). 

XDIF(I)=X(I)-X(I-1),  YD IF(  J  )= Y ( J )- Y( J- 1 ) , 
ZDIF(K)=Z(K)-Z(K- 1) 


Dimensions  in  the  x,  y  and  z  directions  respectively 
associated  with  control  volume  around  the  point 
(X(I),Y(J),Z(K)).  eg.  XCV(I)=XU(I+1)-XU(I).  Thus 
XCV(l)  and  XCV(L1)=0. 
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(O.YL.O) 


(XL.YL.O) 


Fig.  A-1 .  The  computational  domain.  The  coordinate  system  is  right-handed 
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V(I.J.K) 


Node  X(I),Y(J),Z(K) 
is  at  the  geometric  center  of  the 
control  volume;  all  variables 
except  the  velocities  are  stored 
here. 


Fig.  A-2.  The  control  volume. 


APPENDIX  B  -  Program  Source  Code.  The  subroutine  USE  is  for  the 
problem  described  in  Appendix  C. 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C  A  GENERAL  PURPOSE  FORTRAN  PROGRAM  FOR  SOLVING  THREE- 

C  DIMENSIONAL  HEAT  TRANSFER  AND  FLUID  FLOW  IN  RECTANGULAR 

C  COORDINATES 

C 

C  VERSION  1.00,  JULY  1990 

C 

C  DEVELOPED  BY: 

C  S.  B.  SATHE 

C  CODE  69ST ,  DEPARTMENT  OF  MECHANICAL  ENGINEERING 

C  U.S.  NAVAL  POSTGRADUATE  SCHOOL 

C  MONTEREY,  CA  93943 

C  TEL:  ( 408)-646-2417 

C  BITNET  ADDRESS:  5140P@NAVPGS 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

Qk****H********A*****************H***A*********H**H*A* 

LOGICAL  LSTOP 
COMMON/CNTL/LSTOP 

CALL  MESH 
CALL  GEOMET 
CALL  BEGIN 
10  CALL  VARRHO 
CALL  BNDRY 
CALL  PRTOUT 
IF( LSTOP)  STOP 
CALL  COEFF 
GO  TO  10 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  PROFIL 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

COMMON/COEF/FLOW, DIFF, ACOF 
ACOF=DI FF 

IF( FLOW. EQ. 0 . )  RETURN 
TEMP=DI FF-ABS ( FLOW) *0 . 1 
ACOF=Q . 

I F ( TEMP . LE . 0 . )  RETURN 

TEMP-TEMP/DIFF 

ACOF*DI FF*TEMP*  *  5 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  TDMA 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
LOGICAL  LSOLVE,LPRINT,LBLK, LSTOP 

COMMON  F(17,17,13,5),P(17,17,13) ,RHO( 17,17,13) , GAM ( 1 7 , 1 7 , 1 3 ) , 

1  CON( 17,17,13) ,AKP(17, 17, 13) ,AKM( 17, 17 , 13  )  ,AP(17,17,13), 

2  AI P  (  17 , 17 , 1 3  ) ,AIM( 17 , 17 ,13) ,AJP( 17 , 17 , 13) ,AJM( 17,17,13) 
COMMON 

3  X (  1 7 ) , XU (17),XDIF(17) ,XCV( 17 ) ,XCVS( 17 )  , 

4  Y  (  1 7 ) , YV (17) , YDI F ( 17 ) , YCV( 17 ) , YCVS ( 17 ) , 

5  Z  ( 17 ) ,ZW( 17 ) , ZDI F( 17 ) , ZCV( 17 ) ,ZCVS{ 17  )  , 

6  YCVR ( 17 )  , YCVRS{ 17  )  , ARX( 17  )  , ARXJ( 1 7  )  , ARXJP( 17  )  , 

7  R( 17 ) , RMN( 17 ) ,SX( 17 ) , SXMN( 17 ) ,XCVI ( 17 ) ,XCVIP( 17 ) , 

8  YCVJ ( 17 ) , YCVJ P { 17 )  ,ZCVK( 17 ) ,ZCVKP( 17 ) 
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</>•(/>-(/> 


COMMON  DU (17,17,13), DV (17,17,13), DW (17,17,13) , FV (17) ,FVP(17) , 

1  FX( 17 ) , FXM ( 17 ) , FY( 17 ) , FYM( 17 ) , PT { 1 7 ) , QT ( 1 7 ) , 

2  FZ( 17 ) , FZM ( 17 ) 

COMMON/ 1 NDX/RE LAX ( 1 3 ) , LPRINT ( 1 3 )  , LBLK ( 1 1  ) ,NTIMES( 10  )  , 

1LS0LVE( 10) , TIME, DT, XL, YL,ZL, S,RHOCON, ZERO, 

2NF , NFMAX , NP , NRHO , NGAM ,L1,L2,L3,M1,M2,M3,N1,N2,N3, 

3 1 ST, JST, KST, ITER, LAST, 

4IPREF, JPREF, KPREF,MODE 

DIMENSION  D( 17 ) , VAR( 17 ) , VARM( 17 ) , VARP( 17 ) , PHIBAR (17) 

COMMON/HEADIN/TITLE 

CHARACTER* 1 0  TITLE(13) 

QAAAAAAAAAAAAAAAAAAAAAAA'kAAAAAAAA'kAAAAAAAAAAAAAAAAAAAAAA-kA 

I STF= I  ST-1 
JSTF= JST- 1 
KSTF=KST- 1 
IT1=L2+IST 
IT2  =  L3+ I  ST 
JT1=M2+ JST 
JT2=M3+ JST 
KT1=N2+KST 
KT2=N3+KST 


NTSSS=NTIMES(NF) 

DO  999  NT=1 , NTSSS 
DO  391  N=NF , NF 
C 

I F ( . NOT . LBLK ( NF ) ) GO  TO  60 
60  CONTINUE 

COMMENCE  TDMA  LINE  BY  LINE  SWEEPS  FOR  SOLUTION 
DO  90  K=KST,N2 
DO  90  J=JST,M2 
PT ( ISTF ) =0  . 

QT ( I STF ) =F ( ISTF, J,K,N) 

DO  70  1  =  1  ST , L2 

DENOM=AP ( I , J , K ) - PT ( I - 1 ) *AIM( I , J , K ) 

PT ( I ) =AI P ( I , J , K ) /DENOM 

TEMP  =  CON (  I , J , K ) +AJ  P ( I , J , K) *F( I ,  J  +  1  ,  K  ,  N  ) +A JM ( I,J,K)*F(I,J-1,K,N) 
1  +AKP ( I , J,K)*F( I , J , K+ 1 , N ) +AKM ( I,J,K)*F(I,J,K-1,N) 

QT ( I ) = ( TEMP+AIM ( I , J, K) *QT( 1-1 ) )/DENOM 
70  CONTINUE 

DO  80  1 1  =  I  ST , L2 
I-ITl-II 

80  F( I,J,K,N)=F( I  +  1,J,K,N)*PT( I ) +  QT(  I ) 

90  CONTINUE 


DO  190  KK=KST,N3 
K=KT2-KK 

DO  190  JJ= JST , M3 

J»JT2-JJ 

PT ( I STF ) =0  . 

QT( I STF ) =  F ( ISTF, J, K,N) 

DO  170  I  =  I  ST , L2 

DENOM* AP ( I,J,K)-PT(I-1)*AIM(I,J,K) 

PT ( I ) =AI P ( I , J , K ) /DENOM 

TEMP  =  CON (  I , J , K ) +AJP ( I , J , K ) *  F ( I , J+ 1 , K , N ) +AJM (  I,J,K)*F(I,J-1,K,N) 
1  +AKP {  I,J,K)*F( I , J , K  + 1 , N ) +AKM ( I,J,K)*F(I,J,K-1,N) 

QT( I )=( TEMP+AIM ( I , J , K ) *QT ( I - 1 ) ) /DENOM 
170  CONTINUE 

DO  180  1 1  =  I  ST , L2 
I =1 Tl- I I 


180  F(I,J,K,N)=F(I+l,JfK,N)*PT( I ) +QT ( I ) 
190  CONTINUE 


DO  290  I  =  I  ST , L2 
DO  290  K=KST,N2 
PT( JSTF ) =0  . 

QT( JSTF ) =F ( I , JSTF, K,N) 

DO  270  J=JST,M2 

DENOM=AP( I , J , K ) - PT ( J- 1 ) *AJM( I , J, K) 

PT ( J ) =AJP ( I , J,K)/DENOM 

TEMP=CON( I , J , K ) +AKP ( I , J , K) *F( I ,  J  , K+l ,N)+AKM( I/J,K)*F(I,J,K-1,N) 
1  +AI P ( I , J , K ) *  F ( I  +  1 , J , K , N ) + AI M ( I,J,K)*F(I-1,J,K,N) 

QT ( J ) = ( TEMP+AJM ( I , J , K ) *QT( J-l ) ) /DENOM 
270  CONTINUE 

DO  280  J J= JST , M2 
J= JTl— J J 

280  F( I ,  J  , K,N)=F( I ,  J  + 1 , K,N) *  PT  ( J ) +QT ( J  ) 

290  CONTINUE 


DO  390  1 1  =  I  ST , L3 
I=IT2-I I 

DO  390  KK=KST, N3 

K=KT2-KK 

PT( JSTF ) =0  . 

QT ( JSTF ) =F ( I , JSTF , K , N ) 

DO  370  J=JST,M2 

DENOM=AP( I , J, K) -PT( J-l ) *AJM{ I , J, K ) 

PT ( J ) =A  J  P  (  I  ,  J  ,  K  ) /DENOM 

TEMP=CON( I , J , K ) +AKP ( I , J , K ) *  F ( I , J , K  + 1 , N ) +AKM ( I,J,K)*F(I,J,K-1,N) 
1  +AIP{ I , J,K) *F( 1+1 , J,K,N)+AIM( I,J/K)*F(I-1,J,K,N) 

QT ( J ) = ( TEMP+AJM ( I , J , K ) *QT( J-l ) ) /DENOM 
370  CONTINUE 

DO  380  J J=JST , M2 
J=>JTl-JJ 

380  F(I,J,K,N)=F( I,J+1,K,N)*PT(J) +QT( J ) 

390  CONTINUE 


DO  490  J=JST,M2 
DO  490  I  =  1  ST , L2 
PT ( KSTF ) =0 . 

QT( KSTF ) =F ( I , J , KSTF , N ) 

DO  470  K=KST,N2 

DENOM=AP ( I , J , K ) -PT ( K-l ) *AKM( I , J, K) 

PT  (  K  )  =*AKP  (  I  ,  J  ,  K  )  /DENOM 

TEMP-CON( I, J,K)+AIP( I, J,K)*F( 1+1, J,K,N)+AIM( I, J,K)*F( 1-1, J,K,N) 
1  +AJ  P ( I , J, K) *F( I , J+ 1 , K , N ) +AJM ( I,J,K)*F(I,J-1,K,N) 

QT( K ) = ( TEMP+AKM ( I , J , K ) *QT( K-l ) ) /DENOM 
470  CONTINUE 

DO  480  KK=KST,N2 
K=KTl-KK 

480  F( I,J,K,N)=F( I , J , K+l , N ) *PT(K) +QT( K ) 

490  CONTINUE 


DO  590  JJ=JST,M3 
J=JT2-JJ 

DO  590  1 1  =  I  ST , L3 
1= IT2-I I 
PT ( KSTF ) =0  . 

QT ( KSTF ) =F ( I , J , KSTF , N ) 
DO  570  K=KST,N2 
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u  u 


DENOM=AP ( I , J , K ) — PT ( K-l ) * AKH ( I , J, K ) 

PT ( K ) =AKP ( I , J , K ) /DENOM 

TEMP=CON{ I , J , K ) +AI P { I , J, K) *F( I + 1 , J , K , N ) +AIM ( I,J,K)*F(I-1,J,K,N) 
1  + AJ  P ( I,J,K)*F(I,J+1, K , N ) +AJM  ( I,J,K)*F(I,J-1,K,N) 

QT( K ) = ( TEMP+AKM( I , J , K) *QT( K-l ) ) /DENOM 
570  CONTINUE 

DO  580  KK=KST,N2 
K=KTl-KK 

580  F(I,J,K,N)=F(I,J, K+l , N) *PT{ K)+QT( K ) 

590  CONTINUE 

C - 

391  CONTINUE 

999  CONTINUE 

A********************************************** 

ENTRY  RESET 

£*******★**********************★*********★★***** 

DO  400  K=2,N2 
DO  400  J=2 , M2 
DO  400  1=2 , L2 
CON ( I  ,  J , K ) =0 
AP ( I, J,K)=0. 

400  CONTINUE 
RETURN 
END 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  FORM 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
LOGICAL  LSOLVE , L PR I NT , LBLK , LSTOP 

COMMON  F( 17,17,13,5) ,P( 17,17,13) , RHO( 17,17,13) , GAM { 17,17,13), 

1  CON (17, 17, 13) , AKP( 17,17,13) , AKM ( 17 , 1 7 , 1 3 ) ,AP(17,17,13), 

2  AIP(17,17,13) ,AIM( 17,17, 13) ,AJP( 17, 17,13 ) ,AJM( 17,17,13) 

COMMON 

3  X ( 1 7 ) , XU (17) , XDI F ( 1 7 ) ,XCV( 17 ) ,XCVS( 17 ) , 

4  Y ( 1 7  )  , YV (17), YD I F ( 1 7  )  , YCV( 17 )  , YCVS( 17  )  , 

5  Z  ( 17) ,ZW( 17) , ZD I F ( 17) , ZCV( 17 ) , ZCVS( 17 ) , 

6  YCVR ( 17 ) , YCVRS( 17 )  ,ARX( 17 ) , ARX J ( 17 ) ,ARXJP( 17  )  , 

7  R( 17 )  , RMN ( 17 )  , SX( 17 )  , SXMN( 17 )  , XCVI ( 17 )  ,XCVIP( 17  )  , 

8  YCVJ ( 17 ) , YCVJP( 17 ) , ZCVK ( 17 ) , ZCVKP( 17 ) 

COMMON  DU (17, 17, 13), DV (17, 17, 13), DW (17, 17, 13), FV (17), FVP(  17 )  , 

1  FX( 17 ) , FXM( 17 ) , FY( 17 ) , FYM ( 1 7 ) , PT ( 1 7 ) , QT ( 1 7 ) , 

2  FZ ( 1 7 ) , FZM (17) 

COMMON/I NDX/RELAX( 13) ,LPRINT( 13) ,LBLK( 11 ) ,NTIMES( 10)  , 

1 LSOLVE ( 10 ) , TIME , DT ,XL,YL,ZL,S, RHOCON , ZERO , 

2NF , NFMAX , NP , NRHO , NGAM ,Ll,L2,L3,Ml ,M2,M3,Nl,N2,N3, 

3  I  ST , JST , KST , ITER, LAST, 

4IPREF,JPREF,KPREF, MODE 

COMMON/HEADI N/TITLE 
CHARACTER* 1 0  TITLE(13) 

COMMON/CNTL/ LSTOP 
COMMON/SORC/SMAX , SSUM 
COMMON/COEF/FLOW, DIFF, ACOF 

DIMENSION  II (17, 17, 13) ,V(17,17,13) ,W( 17, 17, 13) ,PC( 17,17, 13) 
DIMENSION  T( 17 , 17 , 13  ) 

EQUIVALENCES F(1,1,1,1),U(1,1,1)),(F(1,1,1,2),V(1,1,1)), 

1  ( F( 1  ,1,  1,  3)  ,W( 1, 1, 1 )  )  , <  F(  1  ,  1  ,  1  ,  4  )  , PC( 1  ,  1  ,  1  )  ) 

2  ,(F(1,1,1,5),T(1,1,1)) 

DIMENSION  COFl ( 17 , 17 , 13 , 7 )  , COF 2 ( 1 7 , 1 7  ,  1 3 , 7  )  , COF 3 ( 1 7 , 1 7 , 1 3 , 7  ) 


DIMENSION  COF4 ( 17 , 17 , 1 3 , 7  ) ,  CONl (17,17,13) , CON2 ( 17 , 17 , 1 3 ) 
DIMENSION  CON3( 17,17,13) 

1  FORMAT( 15X, ' COMPUTATION  IN  CARTESIAN  COORDINATES') 

2  FORMAT( 15X, ' COMPUTATION  FOR  AXISYMMETRIC  SITUATION') 

3  FORMAT! 15X, ' COMPUTATION  IN  POLAR  COORDINATES') 

4  FORMAT (14X,40(lH*) ,//) 

C 

COME  HERE  TO  CALCULATE  GRIDS  SPECIFICATION 

eft******************************************************** 
ENTRY  GEOMET 

C 

L2-L1-1 
L3=L2-1 
M2=Ml - 1 
M3=M2-1 
N2=Nl- 1 
N3=N2-1 
X( 1 )=XU( 2  ) 

DO  5  1=2, L2 

5  X(  I ) =  0 . 5 * ( XU ( I  +  1)+XU( I ) ) 

X( LI )=XU( LI ) 

Y ( 1 ) = YV ( 2 ) 

DO  10  J=2 , M2 

10  Y  (  J ) =0 . 5  * ( YV ( J+ 1 ) +  YV ( J )  ) 

Y( Ml ) =YV( Ml  ) 

Z ( 1 ) =ZW ( 2 ) 

DO  7  K=2,N2 

7  Z  (  K  )=0 . 5*  (  ZW(  K+-1  ) +ZW(  K)  ) 

Z(Nl ) =ZW ( Nl ) 


DO  15  1=2, Ll 
15  XDI F ( I ) =X ( I )-X( 1-1  ) 

DO  18  1=2, L2 
18  XCV ( I ) =XU ( I + 1 ) -XU ( I ) 

DO  20  1=3, L2 
20  XCVS( I )=XDIF( I  ) 

XCVS( 3 ) =XCVS ( 3 ) +XD I F ( 2  ) 
XCVS ( L2 ) =XCVS ( L2 ) +XDI F ( Ll ) 
DO  22  1=3, L3 
XCVI ( I ) =  0 . 5*XCV( I ) 

22  XCVI P(I)=XCVI(I) 

XCVIP( 2 )=XCV( 2 ) 

XCVI ( L2 ) =XCV ( L2 ) 

DO  175  K=2 , Nl 
175  ZDIF( K)»Z( K)-Z(K-l) 

DO  178  K  =  2  ,  N2 
178  ZCV( K)=ZW( K+l )-ZW( K) 

DO  270  K=3,N2 
270  ZCVS (K)=ZDIF(K) 

ZCVS( 3 ) =ZCVS ( 3 ) +ZDI F ( 2 ) 
ZCVS( N2 )=ZCVS(N2 )+ZDIF( Nl ) 
DO  272  K=3,N3 
ZCVK ( K ) =0 . 5*ZCV( K ) 

272  ZCVKP( K)=ZCVK( K) 

ZCVKP ( 2 ) =ZCV ( 2 ) 

ZCVK ( N2 ) =ZCV( N2 ) 

DO  35  J  =  2 , M 1 
3  5  Y  D I F ( J ) =  Y ( J ) - Y ( J  - 1  ) 

DO  40  J=2 , M2 


40  YCV ( J ) =YV ( J+ 1 ) - YV ( J ) 

DO  45  J=  3 , M2 
45  Y C V S(J)=YDIF(J) 

YCVS ( 3 ) =YCVS ( 3 ) + YDI F ( 2 ) 

YCVS ( M2 ) =YCVS (M2)+YDIF(Ml ) 

DO  277  J=3 , M3 
YCVJ ( J ) =0 . 5*YCV( J ) 

277  YCVJP ( J ) =YCVJ ( J ) 

YCVJP ( 2 ) =YCV ( 2 ) 

YCVJ ( M2 ) =YCV ( M2 ) 

I F ( MODE . NE . 1 )  GO  TO  55 
DO  52  J=1 , Ml 
RMN (  J  )  =  1 . 0 
52  R ( J ) =1 . 0 
GO  TO  56 

55  DO  50  J=2 , Ml 

50  R ( J ) =R ( J-l )+YDIF(J) 

RMN ( 2 ) =R ( 1  ) 

DO  60  J=3 , M2 

60  RMN( J )=RMN( J-l ) +YCV( J-l ) 

RMN( Ml )=R( Ml ) 

56  CONTINUE 

DO  57  J=1 , Ml 
SX( J)=l. 

SXMN ( J ) =1 . 

I F ( MODE . NE . 3 )  GO  TO  57 
SX ( J ) =R (  J ) 

IF(J.NE.l)  SXMN ( J ) =  RMN ( J ) 

57  CONTINUE 

DO  62  J=2,M2 

YCVR ( J ) =R ( J ) * YCV ( J ) 

ARX ( J ) =  YCVR ( J ) 

I F ( MODE . NE . 3 )  GO  TO  62 
ARX ( J ) =  YCV ( J ) 

62  CONTINUE 

DO  64  J=4,M3 

64  YCVRS( J)=0.5*(R( J)+R( J— 1 ) ) * YDI F ( J ) 
YCVRS ( 3 ) *0 . 5* ( R( 3 ) +R( 1 ) ) *YCVS( 3 ) 
YCVRS( M2 )=0.5*(R(Ml)+R(M3) ) * YCVS (M2 ) 
I F ( MODE . NE . 2  )  GO  TO  67 

DO  65  J=3 , M3 

ARXJ ( J ) =0 . 25* ( 1 . +RMN ( J ) /R ( J ) ) * ARX ( J ) 

65  ARXJP( J ) =ARX( J ) -ARXJ ( J ) 

GO  TO  68 

67  DO  66  J=3,M3 

ARXJ ( J ) =0 . 5*ARX( J ) 

66  ARXJP ( J ) =ARXJ ( J ) 

68  ARXJP ( 2 ) =ARX ( 2 ) 

ARXJ ( M2 ) =ARX ( M2 ) 

DO  70  J=3,M3 

FV ( J ) = ARX J P ( J ) /ARX ( J ) 

70  FVP (J)=1.-FV(J) 

DO  85  1=3, L2 

FX ( I ) =0 . 5*XCV( 1-1 )/XDIF( I ) 

85  FXM ( I )=1 . - FX ( I ) 

FX ( 2 ) =0 . 

FXM { 2  )  =1  . 

FX { LI  )  =1  . 

FXM( Ll ) =0 . 

DO  90  J=  3 , M2 


FY( J ) =0 . 5* YCV( J-l )/YDIF ( 3 ) 

90  FYM{ J ) =1 . -FY{ J ) 

FY ( 2 ) =0  . 

FYM ( 2  )  =1 . 

FY( Ml  )  =1  . 

FYM ( Ml ) =  0  . 

DO  87  K=3,N2 

FZ  (  K ) =0 . 5*ZCV( K-l )/ZDIF( K) 

87  FZM(K)=1.-FZ(K) 

FZ ( 2  )  =0  . 

FZM ( 2 ) =1 . 

FZ ( Nl ) =1  . 

FZM(Nl ) =  0  . 

CON , AP , U , V , RHO , PC  AND  P  ARRAYS  ARE  INITIALIZED  HERE 
DO  95  K=1 , Nl 
DO  95  J  =  1  , Ml 
DO  95  I =1 , Ll 
PC( I , J , K  )  =0  . 

U( I , J,K)=0. 

V  (  I  ,  J , K ) =0 . 

W( I , J , K  )  =0  . 

CON ( I , J,K)=0. 

AP ( I , J,K)=0. 

RHO ( I , J , K ) =RHOCON 
P ( I , J , K ) =0  . 

95  CONTINUE 

I F ( MODE . EQ . 1 )  PRINT  1 
I F ( MODE . EQ . 2 )  PRINT  2 
I F ( MODE . EQ . 3 )  PRINT  3 
PRINT  4 
RETURN 


C 

COME  HERE  TO  CALCUALTE  COEFFICIENTS  FOR  FINITE  DIFF.  EQNS . 


:*****************] 


******** 


ENTRY  COEFF 


£******************************* 


*************** 


c 

COEFFICIENTS  FOR  THE  U  EQUATION 
C 


CALL  RESET 
NF=1 

I F (  . NOT . LSOLVE ( NF ) )  GO  TO  100 

I  ST=  3 

JST=2 

KST=2 

CALL  DIFFUS 
REL=1 . -RELAX ( NF ) 

DO  103  1=3, L2 
DO  103  J=2 , M2 
DO  103  K=2 , N2 

COEFFICIENTS  AEAST  AND  AWEST 

FL=U ( I , J , K ) * { FX ( I ) *RHO( I , J , K ) +FXM( I ) *  RHO ( I - 1 , J , K )  ) 

FLP  =  U ( I+1,J,K)*(FX(I+1) *  RHO ( I  + 1 , J , K )  +  FXM (  1  +  1 ) *RHO( I , J , K)  ) 

FLOW* YCV { J ) *ZCV ( K ) *  0 . 5  * ( FL+  FLP ) 

DIFF=YCV( J ) *ZCV ( K ) *GAM ( I , J , K)/XCV( I ) 

CALL  PROFIL 

AIM( 1+1 , J,K) =ACOF+AMAXl( ZERO, FLOW) 

AIP( I , J,K)=AIM( 1+1, J,K)-FLOW 
COEFFICIENTS  ANORTH  AND  ASOUTH 

FL-XCVI ( I ) * V ( I , J  + 1 , K ) * ( F Y ( J  +  1 ) *  RHO ( I , J  +  l , K ) +  FYM ( J  +  1 ) *RHO(  I , J , K )  ) 
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FLM=XCVI P ( 1-1 ) *V{ 1-1 ,J  +  1 , K ) *  (  FY{  J+l ) *RHO(  1-1 ,  J  +  l , K ) +  FYM(  J  +  l ) * 

1  RHO ( 1-1, J,K) ) 

GM=GAM( I , J , K ) *GAM( I , J+l  ,  K) 

1  /  (  YCV  (  J  )  *GAM  (  I  ,  J  +  l  ,  K-)+YCV(  J  +  l  )  *GAM(  I  ,  J  ,  K  )  + 

2  1 . 0E-30 ) *XCVI ( I  ) 

GMM»GAM{ 1-1 , J, K) *GAM( 1-1 , J+l , K) 

1  / ( YCV ( J ) *GAM( 1-1 , J+l , K ) + YCV ( J+l ) * 

2  GAM( I-l,J,K)+l.E-30 ) *XCVIP( 1-1 ) 

DIFF=ZCV( K) *2 . * ( GM+GMM ) 

FLOW=ZCV( K) * ( FL+FLM) 

CALL  PROFIL 

AJM (  I  ,  J  +  l  ,  K  )  =ACOF+AMAXl  (  ZERO,  FI.OW) 

AJP  ( I , J , K)=AJM( I , J+l , K)-FLOW 
COEFFICIENTS  AIN  AND  AOUT 

FL=XCVI ( I ) *W( I , J , K+l ) * ( FZ( K+l ) *RHO( I , J , K+l )+FZM( K+l ) *RHO( I , J , K ) ) 
FLM=XCVIP( I-1)*W(I-1,J,K  +  1)*(FZ(K  +  1) *RHO( 1-1 , J , K  +  l )+FZM( K  +  l  )  * 

1  RHO ( 1-1, J,K) ) 

GM=GAM( I , J , K ) *GAM{ I , J , K+l ) 

1  / ( ZCV ( K ) *GAM { I , J , K+l )+ZCV( K+l ) *GAM ( I , J , K ) + 

2  1 . 0E-30 ) *XCVI ( I ) 

GMM=GAM( 1-1 , J, K) *GAM( 1-1 , J, K+l ) 

1  / ( ZCV ( K ) *GAM ( 1-1 , J , K+l ) +ZCV( K+l ) * 

2  GAM ( 1-1 ,J, K ) + 1 . E-30 ) *XCVIP( 1-1 ) 

DIFF=YCV( J ) *2 . * ( GM+GMM) 

FLOW=YCV ( J ) * ( FL+FLM) 

CALL  PROFIL 

AKM ( I , J , K+l ) =ACOF+AMAXl ( ZERO , FLOW ) 

AK  P ( I, J,K) =AKM ( I , J , K  +  l )-FLOW 
103  CONTINUE 

DO  104  J=2 , M2 
DO  104  K=2 , N2 

COEFFICIENTS  AEAST  AND  AWEST 
ARE A= YCV ( J ) *  ZCV ( K ) 

FLOW=AREA*RHO( 1,J,K)*U(2,J,K) 

DIFF=AREA*GAM( 1 , J,K)/XCV( 2) 

CALL  PROFIL 

AIM ( 3 , J , K ) =ACOF+ AMAXl ( ZERO , FLOW ) 

FLOW=AREA*RHO( Ll ,J,K)*U(Ll,J,K) 

DI FF=AREA*GAM( Ll , J , K )/XCV ( L2 ) 

CALL  PROFIL 

AIP ( L2 , J , K ) =ACOF+ AMAXl ( ZERO , FLOW ) -FLOW 

104  CONTINUE 

DO  105  1=3, L2 
DO  105  K=2,N2 

COEFFICIENTS  ANORTH  AND  ASOUTH 

FL=XCVI ( I ) * V ( I , 2,K) *RHO( I , 1 , K) 

FLM=XCVI P ( 1-1 ) *V( 1-1 , 2, K) *RHO( I - 1 , 1 , K  ) 

FLOW=ZCV ( K ) * ( FL+FLM) 

GM=XCVI ( I ) *GAM ( I , 1 , K ) +XCVIP( 1-1 ) *GAM ( 1-1 , 1 , K ) 
DIFF=ZCV(K)*GM/YDIF( 2) 

CALL  PROFIL 

AJM ( I , 2 , K ) =ACOF+ AMAXl (ZERO, FLOW) 

FL=XCVI ( I ) * V ( I , Ml , K ) *RHO ( I , Ml , K ) 

FLM=XCVI P ( I-l)*V(I-l,Ml,K) *RHO( 1-1 , Ml , K ) 

FLOW=ZCV ( K ) * ( FL+FLM) 

GM=XCV I ( I ) *GAM ( I , Ml , K)+XCVIP( 1-1 ) *GAM ( I - 1 , Ml , K ) 

DI FF=ZCV ( K) *GM/YDI F ( Ml ) 

CALL  PROFIL 

AJP( I , M2 , K ) =ACOF+ AMAXl ( ZERO, FLOW) -FLOW 

105  CONTINUE 
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DO  106  1=3, L2 
DO  106  J=2 , M2 
COEFFICIENTS  AIN  AND  AOUT 

FL=XCVI (  I ) *W( I , J , 2 ) *  RHO ( I , J , 1  ) 

FLM=XCV I P ( 1-1 ) *W( 1-1 , J , 2 ) *RHO( 1-1 , J , 1 ) 

FLOW=YCV ( J ) * ( FL+FLM) 

GM=XCVI ( I ) *GAM(  I , J , 1 ) +XCVIP( 1-1 ) *  GAM ( I - 1 , J  ,  1  ) 

DIFF=YCV( J ) *GM/ZDIF( 2 ) 

CALL  PROFIL 

AKM ( I , J , 2 ) =ACOF+AMAXl ( ZERO, FLOW) 

FL=XCVI ( I ) *W( I , J , Nl ) *RHO ( I , J , Nl  ) 

FLM=XCVIP( 1-1 ) *W( 1-1 ,J,Nl ) *  RHO { 1-1 , J,Nl  ) 

FLOW=YCV ( J ) * ( FL+FLM) 

GM=XCVI ( I ) *GAM ( I , J,Nl )+XCVIP( 1-1 ) *GAM ( 1-1 , J , Nl  ) 

DIFF=YCV<  J ) *GM/ZDI F ( Nl ) 

CALL  PROFIL 

AKP ( I , J , N2 ) =ACOF +AMAX1 ( Z ERO , FLOW )- FLOW 
106  CONTINUE 

DO  107  1=3, L2 
DO  107  J=2 , M2 
DO  107  K=2,N2 
VOL= YCV ( J ) *XCVS( I ) *ZCV ( K ) 

APT= ( RHO ( I , J , K ) *XCVI ( I ) +RHO( I - 1 , J , K ) *XCVI P  (  1-1 )  ) 

1/(XCVS( I ) *DT ) 

AP ( I , J , K ) =AP ( I , J , K ) -APT 

CON ( I , J, K)=CON( I , J , K ) +APT*U ( I , J , K ) 

AP ( I , J , K ) = 

1  ( -AP ( I , J , K ) *VOL+AIP( I , J, K)+AIM{ I , J , K ) +A J  P ( I , J , K ) +AJM ( I , J, K ) 

2  +AKP( I , J, K)+AKM( I , J, K) ) 

3/RELAX (NF) 

CON ( I , J , K ) =CON ( I , J , K ) * VOL+REL*AP ( I,J,K)*U(I,J,K) 

DU ( I , J , K ) =VOL/XDI F  (  I ) 

DU ( I , J, K)=DU( I , J,K)/AP( I , J, K) 

107  CONTINUE 

DO  7099  1=1, Ll 
DO  7099  J=1 , Ml 
DO  7099  K=1 , Nl 
COFl ( I , J , K, 1 )=AIP( I , J, K ) 

COFl ( I , J,K,2)=AIM( I , J , K ) 

COFl ( I , J , K , 3 ) =AJP ( I , J , K) 

COFl ( I , J , K, 4 ) =AJM ( I , J, K) 

COFl ( I , J , K, 5 )=AKP( I , J, K) 

COFl ( I , J , K , 6 ) =AKM( I,J,K) 

COFl(I,J,K,7) =AP ( I , J , K ) 

CONl( I , J , K ) =CON ( I , J , K  ) 

7099  CONTINUE 

C  TEMPORARY  USE  OF  PC(I,J)  TO  STORE  UHAT 

DO  151  K=2,N2 
DO  151  J=2 , M2 
DO  151  1=3, L2 

151  PC(I,J,K)=(AIP(I,J,K)*U(I+1,J,K) +AI M ( I,J,K) *U( I - 1 , J , K ) 

1  +AJP ( I , J, K) *U( I , J+l , K) +AJM( I , J , K ) *U( I , J- 1 , K ) 

2  +AKP ( I , J,K)*U( I , J , K+ 1 ) +AKM ( I,J,K)*U(I,J,K-1) 

3  +CON( I, J,K) )/AP( I, J,K) 

100  CONTINUE 

C 

COEFFICIENTS  FOR  THE  V  EQUATION - 

C 

CALL  RESET 
NF  =  2 
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IF( . NOT . LSOLVE ( NF ) )  GO  TO  200 
I  ST=2 
JST=  3 
KST=2 

CALL  DIFFUS 
REL=1 . -RELAX ( NF ) 

DO  203  J=3 , M2 
DO  203  1=2, L2 
DO  203  K=2 , N2 

COEFFICIENTS  ANORTH  AND  ASOUTH 

FL= V  (  I , J , K ) * ( F Y ( J ) *  RHO ( I , J , K ) +  F YM ( J  )  *  RHO  (  I  ,  J - 1 ,  K  )  ) 

FLP-V ( I , J  +  l , K) *( FY{ J  +  l ) *  RHO ( I ,  J  +  l , K ) +FYM ( J  +  l ) *RHO( I , J , K )  ) 
FLOW=XCV( I ) *  ZCV { K ) *  0 . 5  * ( FL+FLP ) 

DIFF=XCV( I ) *ZCV ( K ) *GAM ( I , J, K )/YCV( J ) 

CALL  PROFIL 

AJM  (  I , J  + 1 , K ) -ACOF+AMAXl ( ZERO , FLOW ) 

AJP ( I , J , K )=AJM( I , J+l , K)-FLOW 
COEFFICIENTS  AEAST  AND  AWEST 

FL- YCVJ  (  J  )  *U ( I +1 , J , K ) * ( FX( 1+1 )  *RHO(  1  +  1  ,  J  ,  K  ) +FXM(  1  +  1  )  *RHO(  I  ,  J  ,  K)  ) 
FLM- YCVJP ( J-l )  *U  ( 1  +  1 ,J-1 ,K) *( FX( 1  +  1 ) *RHO( 1  +  1 , J-l , K)+FXM(  1  +  1 ) * 

1  RHO ( I , J-l , K ) ) 

GM=GAM( I , J , K ) *GAM ( 1+1 , J,K) 

1  /(XCV( I )*GAM( 1+1, J,K)+XCV< I+1)*GAM( I,J,K)+ 

2  1 . 0E-30 ) * YCVJ ( J ) 

GMM=GAM( I , J-l ,K) *GAM( 1+1 , J-l , K) 

1  /( XCV( I ) *GAM( 1+1 , J-l , K)+XCV( 1+1 ) * 

2  GAM( I,J-l,K)+l.E-30) * YCVJP ( J-l ) 

DI FF=ZCV ( K ) *  2 . * ( GM+GMM ) 

FLOW-ZCV ( K ) * ( FL+FLM) 

CALL  PROFIL 

AIM ( 1+1 , J,K) =ACOF+AMAXl( ZERO, FLOW) 

AIP( I , J , K ) =AIM( 1+1, J,K)-FLOW 
COEFFICIENTS  AIN  AND  AOUT 

FL=YCVJ( J) *W( I , J,K  +  1 ) * ( FZ( K  +  l ) *RHO( I , J , K+ 1 )  +  FZM( K+l ) *RHO(  I , J , K )  ) 
FLM= YCVJP (J-l ) *W( I , J-l , K+l ) * ( FZ ( K+l ) *RHO( I , J-l , K+l )+FZM( K+l ) * 

1  RHO ( I , J-l , K  )  ) 

GM=GAM( I , J , K ) *GAM ( I , J, K+l  ) 

1  /( ZCV{ K ) *GAM( I , J , K+l ) +ZCV( K+l ) *GAM ( I , J , K ) + 

2  1 . 0E-30 ) * YCVJ ( J ) 

GMM=GAM ( I , J-l , K ) *GAM ( I , J-l , K+l ) 

1  /( ZCV ( K ) *GAM( I , J-l , K+l ) +ZCV( K+l ) * 

2  GAM ( I,J-1,K)+1.E-30)*YCVJP( J-l) 

DI FF=XCV ( I ) *2 . * (GM+GMM) 

FLOW=XCV ( I ) * ( FL+FLM) 

CALL  PROFIL 

AKM( I , J, K+l ) =ACOF+AMAXl ( ZERO, FLOW) 

AKP ( I , J , K ) =AKM ( I , J, K+l )-FLOW 
203  CONTINUE 

DO  204  1=2, L2 
DO  204  K=2 , N2 

COEFFICIENTS  ANORTH  AND  ASOUTH 
AREA=XCV( I ) *ZCV ( K ) 

FLOW-AREA* RHO ( I , 1 , K ) * V ( I , 2 , K ) 

DI FF=AREA*GAM ( I , 1 , K)/YCV( 2 ) 

CALL  PROFIL 

AJM ( I , 3 , K ) -ACOF+AMAXl ( ZERO, FLOW) 

FLOW=AREA*RHO( I , Ml , K ) *V( I , Ml , K ) 

D I FF -AREA  * GAM ( I , Ml , K )/YCV( M2 ) 

CALL  PROFIL 

AJP ( I , M2 , K ) -ACOF+AMAXl ( Z  ERO , FLOW ) - FLOW 
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204  CONTINUE 

DO  205  J=3 , M2 
DO  205  K=2,N2 

COEFFICIENTS  AEAST  AND  AWEST 

FL= YCVJ ( J ) *U ( 2 , J , K ) *  RHO ( 1 , J , K ) 

FLM= YCVJP (J-1)*U(2,J-1,K) *RHO( 1,J-1,K) 

FLOW=ZCV( K) * ( FL+FLM) 

GM-YCVJ ( J ) *GAM ( 1 ,  J  ,  K  ) +YCVJP ( J-l  ) *  GAM { 1 , J-l , K ) 

DIFF=ZCV{ K) *GM/XD I F ( 2 ) 

CALL  PROFIL 

AIM( 2 , J , K ) =ACOF+AMAXl ( ZERO , FLOW ) 

FL=YCVJ( J)*U(H,  J,K)*RHO(Ll,  J,K) 

FLM= YCVJP (J— l)*U(Ll,J— 1,K) *RHO( Ll,J— 1,K) 

FLOW—ZCV ( K) * ( FL+FLM) 

GM* YCVJ ( J ) *GAM ( LI , J , K ) + YCVJP ( J-l ) *  GAM ( Ll , J-l , K ) 

DI FF=ZCV ( K ) *GM/XD I F ( LI ) 

CALL  PROFIL 

A I P ( L2 , J , K ) =ACOF  +AMAX1 ( ZERO , FLOW) -FLOW 

205  CONTINUE 

DO  206  1=2, L2 
DO  206  J=  3 , m2 
COEFFICIENTS  AIN  AND  AOUT 

FL  =  YCVJ( J ) *W( I , J, 2 ) *RHO ( I , J, 1  ) 

FLM= YCVJP (J-l ) *W( I , J-l , 2 ) *RHO ( I , J-l  ,  1  ) 

FLOW=XCV ( I ) * ( FL+FLM ) 

GM=YCVJ( J ) *GAM ( I , J , 1 ) + YCVJP ( J-l ) *GAM ( I , J-l , 1 ) 

DI FF=XCV ( I ) +GM/ZDI F ( 2 ) 

CALL  PROFIL 

AKM ( I , J , 2 ) =ACOF +AMAX1 ( ZERO , FLOW ) 

FL=YCVJ( J ) *W( I , J , N1 ) *  RHO ( I , J ,Nl ) 

FLM= YCVJP ( J-l ) *W( I , J-l ,N1 ) *  RHO ( I , J-l ,Nl ) 

FLOW=XCV( I ) * ( FL+FLM) 

GM=YCVJ ( J ) *GAM ( I , J , Nl ) + YCVJP ( J-l ) *GAM{ I , J-l , Nl ) 

DI FF=XCV ( I ) *GM/ZDI F ( Nl ) 

CALL  PROFIL 

AKP ( I , J ,N2 )=ACOF+AMAXl (ZERO, FLOW) -FLOW 

206  CONTINUE 

DO  207  1=2 , L2 
DO  207  J=3 , M2 
DO  207  K=2 , N2 
VOL  =  XCV( I ) *  YCVS ( J ) *ZCV( K ) 

APT= ( RHO (  I , J , K ) *  YCVJ ( J ) +RHO(  I , J-l , K ) *YCVJP( J-l )  ) 

1/ ( YCVS ( J ) *DT ) 

AP ( I , J , K ) «AP ( I , J , K ) -APT 

CON ( I , J , K ) =CON ( I , J,K)+APT*V( I , J , K ) 

AP( I , J,K)= 

1  ( -AP ( I, J,K)*VOL+AIP( I , J , K ) +AIM ( I , J , K ) +AJP ( I , J , K ) +AJM ( I, J,K) 

2  +AKP ( I , J , K ) +AKM ( I , J , K ) ) 

3/RELAX( NF ) 

CON ( I,J,K)«CON( I ,J , K) *VOL+REL*AP( I , J , K ) * V ( I , J , K ) 

DV ( I , J , K ) = VOL/ YD I F ( J ) 

DV( I, J,K)=DV( I , J,K)/AP( I , J,K) 

207  CONTINUE 

DO  0099  1=1, Ll 
DO  8099  J= 1 , Ml 
DO  8099  K= 1 , Nl 
COF2 ( I , J , K , 1 ) =AI P ( I , J , K  ) 

COF2 ( I , J , K , 2 ) =AIM( I , J , K ) 

COF2 ( I , J, K, 3 ) =AJP ( I , J, K) 

COF2 ( I,J,K,4)=AJM( I, J,K) 
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COF2(I,J,K,5)=AKP(I,J,K) 

COF2(I,J,K,6)=AKM(I,J,K) 

C0F2 ( I ,  J  ,  K  ,  7  ) =AP ( I, J,K) 

C0N2 ( I , J , K ) -CON ( I , J , K ) 

8099  CONTINUE 
200  CONTINUE 
C 

COEFFICIENTS  FOR  THE  W  EQUATION - 

C 

CALL  RESET 
NF-3 

IF( . NOT . LSOLVE ( NF ) )  GO  TO  300 
I  ST-2 
JST-2 
KST-  3 

CALL  DIFFUS 
REL-1 . -RELAX ( NF ) 

DO  303  K-3 , N2 
DO  303  J-2 , M2 
DO  303  1=2, L2 

COEFFICIENTS  AIN  AND  AOUT 

FL=W ( I , J , K ) * ( FZ ( K ) *RHO( I , J, K ) +  FZM ( K ) *RHO( I , J , K— 1  )  ) 

FLP=W( I , J , K  +  l ) * ( FZ ( K  +  l ) *RHO( I , J , K+ 1 ) +FZM ( K+ 1 ) *RHO( I , J , K)  ) 
FLOW=YCV ( J ) *XCV (I)*0.5*(FL+FLP) 

DIFF=YCV( J)*XCV( I )*GAM( I , J , K )/ZCV( K ) 

CALL  PROFIL 

AKM ( I , J, K+l )=ACOF+AMAXl ( ZERO, FLOW) 

AKP ( I , J , K ) -AKM ( I , J , K+l ) -FLOW 
COEFFICIENTS  ANORTH  AND  ASOUTH 

FL=ZCVK( K ) *V( I , J+l , K ) * ( FY( J+l ) *RHO< I , J+l , K ) +FYM ( J+ 1 ) *RHO( I , J , K ) ) 
FLM=ZCVKP( K-l ) *V( I , J+l , K-l ) * ( FY( J+l ) *RHO( I , J+l , K-l ) +FYM( J+l ) * 

1  RHO ( I , J, K-l  )  ) 

GM=GAM( I , J , K ) *GAM ( I , J+l ,K) 

1  / ( YCV ( J ) *GAM ( I , J+l , K)+YCV( J+l ) *GAM( I , J , K ) + 

2  1 . 0  E- 3  0 ) *  ZCVK ( K ) 

GMM=GAM( I , J , K-l ) *GAM( I , J+l , K-l ) 

1  / ( YCV ( J ) *GAM ( I , J+l , K-l ) +YCV( J+l ) * 

2  GAM ( I , J , K-l ) +1 . E-30 ) *ZCVKP(  K-l ) 

DIFF=XCV( I ) *2 . * ( GM+GMM ) 

FLOW=XCV( I ) * ( FL+FLM) 

CALL  PROFIL 

AJM ( I , J+l ,K)=ACOF+AMAXl( ZERO, FLOW) 

A J  P ( I , J , K ) -  A  J  M ( I , J+l , K) -FLOW 
COEFFICIENTS  AEAST  AND  AWEST 

FL-ZCVM  K) *U( I +1 , J , K ) * ( FX( 1+1) *RHO( I + 1 , J , K ) + FXM ( 1+1 ) *RHO( I , J , K) ) 
FLM«ZCVKP< K-l ) *U( I + 1 , J , K-l ) *( FX { I + 1 ) * RHO ( 1+1 , J , K-l ) +FXM( 1+1 ) * 

1  RHO( I , J , K-l  )  ) 

GM=GAM ( I , J , K ) *GAM ( I+1,J,K) 

1  /( XCV( I ) *GAM{ 1+1 , J , K ) +XCV{ 1+1 ) *GAM ( I , J , K ) + 

2  1 . 0E- 30 ) *ZCVK ( K ) 

GMM=GAM( I , J, K-l ) *GAM ( 1+1 , J, K-l ) 

1  /( XCV( I ) *GAM ( 1+1 , J , K-l ) +XCV( 1+1 ) * 

2  GAM ( I , J , K-l )+l .E-30) *  ZCVK P (K-l ) 

DIFF=YCV( J ) *2 . * ( GM+GMM) 

FLOW- YCV ( J) M FL+FLM) 

CALL  PROFIL 

AIM ( I + 1 , J , K ) -ACOF+AMAXl ( ZERO, FLOW) 

AIP( I , J , K ) -AIM ( I + 1 , J , K ) -FLOW 
303  CONTINUE 

DO  304  J-2 , M2 
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DO  304  1=2, L2 

COEFFICIENTS  AIN  AND  AOUT 
AREA=YCV<  J ) *XCV( I ) 

FLOW=AREA*RHO ( I  ,  J  ,  1 ) *W( I , J , 2  ) 

DIFF=AREA*GAM( I , J , 1 ) /ZCV ( 2 ) 

CALL  PROFIL 

ARM ( I , J , 3 ) =ACOF+AMAXl ( ZERO , FLOW ) 

FLOW=AREA*  RHO ( I , J , Nl ) *W ( I , J , Nl ) 

DIFF=AREA*GAM( I , J , Nl ) /ZCV ( N2 ) 

CALL  PROFIL 

AKP ( I , J ,N2 )=ACOF+AMAXl ( ZERO, FLOW) -FLOW 

304  CONTINUE 

DO  305  1=2, L2 
DO  305  K=  3 , N2 

COEFFICIENTS  ANORTH  AND  ASOUTH 

FL=ZCVK ( K ) *V( I , 2 , K ) *  RHO ( I , 1 , K) 

FLM=ZCVKP( K-l ) *V( I , 2 , K-l ) *RHO( 1,1, K-l ) 

FLOW=XCV ( I ) * ( FL+FLM ) 

GM=ZCVK( K ) *GAM( I , 1 , K ) +ZCVKP ( K-l ) *GAM( I , 1 , K-l ) 

DIFF=XCV( I ) *GM/YDI F ( 2 ) 

CALL  PROFIL 

AJM( 1,2, K )=ACOF+AMAXl (ZERO, FLOW) 

FL  =  ZCVK ( K ) *  V (  I , Ml , K ) *RHO( I,Ml,K) 

FLM=ZCVKP( K-l ) *V( I , Ml , K-l ) *RHO( I ,Ml , K-l ) 

FLOW=XCV( I ) * ( FL+FLM) 

GM=ZCVK ( K) *GAM ( I ,Ml , K)+ZCVKP( K-l ) *GAM ( I ,Ml , K-l ) 

DIFF=XCV( I ) *GM/YDI F ( Ml ) 

CALL  PROFIL 

AJ  P ( I , M2 , K)=ACOF+AMAXl ( ZERO, FLOW ) -FLOW 

305  CONTINUE 

DO  306  K=3 , N2 
DO  306  J=2 , M2 

COEFFICIENTS  A EAST  AND  AWEST 

FL=ZCVK ( K ) *U( 2 , J , K ) *RHO( 1 , J , K ) 

FLM=ZCVKP( K-l ) *U( 2 , J , K-l ) *RHO( 1 , J , K-l ) 

FLOW= YCV ( J ) * ( FL+FLM) 

GM=ZCVK ( K ) *GAM ( 1 , J, K ) +ZCVKP( K-l ) *GAM ( 1 , J , K-l ) 

DIFF=YCV( J ) *GM/XDI F ( 2 ) 

CALL  PROFIL 

AIM( 2 , J , K ) =ACOF +AMAX1 (ZERO, FLOW) 

FL-ZCVK ( K ) *U ( Ll , J , K ) *RHO( Ll , J , K ) 
FLM=ZCVKP(K-l)*U(Ll,J,K-l)*RHO(Ll,J,K-l) 

FLOW- YCV ( J ) * { FL+FLM ) 

GM-ZCVK ( K ) *GAM ( Ll , J , K ) +ZCVKP (K-l ) *GAM ( Ll , J , K-l ) 

DI FF-YCV ( J ) *GM/XDI F ( Ll ) 

CALL  PROFIL 

AIP(L2,J,K)=ACOF+AMAXl( ZERO, FLOW) -FLOW 

306  CONTINUE 

DO  307  1-2, L2 
DO  307  J-2 , M2 
DO  307  K=  3 , N2 
VOL-YCV ( J ) *  ZCVS ( K ) *XCV( I ) 

APT- ( RHO ( I , J , K ) *ZCVK ( K ) +RHO ( I , J , K- 1 ) *  ZCVKP ( K- 1 )  ) 

1/ ( ZCVS ( K ) *  DT ) 

AP ( I , J , K ) =AP ( I , J , K ) -APT 

CON ( I , J , K ) =CON ( I , J , K ) +APT*W ( I,J,K) 

AP ( I , J , K ) = 

1  ( - AP ( I , J , K ) *VOL+AI P ( I,J,K)+AIM(I,J,K)fAJP(I,J,K) +AJM( I , J , K ) 

2  +AKP( I , J, K )+AKM( I , J , K) ) 

3/RELAX ( NF ) 
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CON ( I , J , K ) =CON ( I , J , K ) * VOL+REL*AP ( I , J , K ) *W (  I  ,  J  ,  K  ) 

DW( I , J , K ) = VOL/ZD I F(  K  ) 

DW ( I , J , K ) =DW ( I , J , K )/AP { I , J , K) 

307  CONTINUE 

DO  9099  1=1 , Ll 
DO  9099  J=1 , Ml 
DO  9099  K=1 , Nl 
COF3(I,J,K,l) =AI P ( I , J , K  ) 

COF 3 ( I , J , K , 2 ) =AIM ( I , J , K ) 

COF  3 ( I , J , K , 3 ) =AJ  P ( I , J, K) 

COF3( I , J, K, 4 )=AJM( I , J , K ) 

COF3 ( I , J , K , 5 ) =AKP ( I , J, K) 

COF3(I,J,K,6) =AKM ( I , J , K ) 

COF  3 ( I, J, K, 7 )=AP( I , J, K) 

CON  3 ( I , J , K ) =CON ( I , J , K ) 

9099  CONTINUE 
300  CONTINUE 
C 

COEFFICIENTS  FOR  THE  PRESSURE  EQUATION _ 

C  ' 

NF=NP 

I F ( . NOT .LSOLVE(NF) )  GO  TO  500 
I  ST=2 
J  ST=2 
KST=2 

DO  501  J=2 , M2 
DO  501  K=2 , N2 
AIM( 2 , J , K  )  =0 . 0 
AIP(L2, J#K)=0.0 

CON (  2  ,  J  ,  K  ) =RHO (1,J,K)*U{2,J,K) *  YCV ( J ) *  ZCV ( K ) 

501  CONTINUE 

DO  502  1=2, L2 
DO  502  K=2 , N2 
AJM ( I , 2 , K  )=0 . 0 
A J P ( I , M2 , K ) =  0 . 0 

CON ( I , 2, K)=RHO{ I , 1 , K) *V( I , 2,K) *XCV( I ) *ZCV ( K ) 

502  CONTINUE 

DO  503  1=2, L2 
DO  503  J=2 , M2 
AKM ( I , J , 2 ) =0 . 0 
AKP ( I , J , N2  )  =0 . 0 

CON ( I , J, 2 )=RHO( I , J , 1 ) *W ( I , J , 2 ) *XCV( I ) *  YCV ( J ) 

503  CONTINUE 
C 

DO  504  K=2 , N2 
DO  504  J=2 , M2 
DO  504  1=2, L2 
AREA-YCV ( J ) *  ZCV ( K ) 

ARHO-AREAM  FX( 1  +  1 ) *RHO{ 1  +  1 , J , K ) +FXM{ 1  +  1 ) *RHO(  I  ,  J,  K)  ) 
FLOW=ARHO*  PC ( 1  +  1, J,K) 

I F ( I . EQ . L2 ) FLOW=ARHO*U ( Ll , J, K) 

AIP( I , J,K)=ARHO*DU( 1+1, J,K) 

AIM( I + 1 , J , K ) =AI P ( I, J,K) 

CON ( I , J , K ) =CON ( I , J , K ) - F  LOW 
CON ( 1  +  1  ,  J  ,  K) =CON( 1  +  1 , J , K ) +  FLOW 
C 

AREA=XCV (  I  ) *  ZCV ( K ) 

ARHO=AREA* ( FY( J+l ) *RHO( I , J  + 1 , K ) +  FYM ( J+ 1 ) * RHO(  I , J , K  )  ) 
VHAT  =  (COF2(  I  ,  J  +  l , K, 1 ) *V( 1  +  1 , J  +  l , K) 

1  +COF2(I,J+l,K,2)*V(I-l,J+l,K) 
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2  +COF2(I,J+l,K,3)*V(I,J+2,K) 

3  +C0F2( I , J+l , K , 4 )*V( I,J,K) 

4  +  C0F2( I , J  +  l , K, 5 ) *V< I , J  +  l , K  +  l  ) 

5  +COF2(I,J+l,K,6)*V{I,J+l,K-l) 

6  +C0N2 ( I , J+l , K) )/( C0F2( I,J+l,K,7)+1.0E-30) 
FLOW=ARHO* VHAT 

I F ( J . EQ . M2 ) FLOW=ARHO  * V ( I,Ml,K) 

AJ P ( I , J, K)=ARHO*DV( I , J  +  l , K  ) 

AJM ( I , J+l ,K)=AJP( I  ,  J  ,  K  ) 

CON ( I , J , K ) =  CON ( I  ,  J  ,  K  ) - FLOW 
CON ( I , J  +  l , K)=CON( I , J  +  l , K ) +  FLOW 

AREA=XCV( I ) * YCV ( J ) 

ARHO=AREA* ( FZ ( K+l ) *RHO( I , J , K+l ) +FZM( K+l ) *RHO( I , J, K)  ) 
WHAT  =  ( COF3 ( I , J, K+l , 1 ) *W( 1+1 , J , K+l ) 

1  +COF3( I , J, K+l , 2) *W( 1-1 , J, K+l ) 

2  +COF3( I , J, K  +  l , 3 ) *W( I , J  +  l  ,  K+l ) 

3  +COF3 ( I , J , K+l , 4 ) *W( I , J-l , K+l ) 

4  +COF3( I , J, K+l , 5 ) *W( I ,  J  ,  K  +  2 ) 

5  +COF3( I,J,K+1,6)*W( I, J,K) 

6  +CON3 ( I , J , K+ 1 ) )/{ COF3 ( I,J,K+l,7)+l.E-30) 
FLOW=ARHO*WHAT 

I F ( K. EQ.N2 ) FLOW=ARHO*W( I , J , Nl ) 

AKP ( I , J , K ) = ARHO  *  DW ( I , J , K  +  l ) 

AKM ( I , J, K  +  l )=AKP( I  ,  J  ,  K  ) 

CON ( I ,  J,K)=CON( I , J,K)-FLOW 
CON ( I , J , K+l  ) =FLOW 

AP ( I , J , K  )  =  AI P ( I , J , K ) +AIM ( I , J , K ) +AJ P ( I , J , K ) +AJM ( I , J  ,  K ) 

1  +AKP(I,J,K)+AKM(I,J,K) 

504  CONTINUE 

DO  703  1=1, LI 
DO  703  J= 1 , Ml 
DO  703  K=1 , Nl 
COF  4 ( I , J , K , 1 ) =AI P ( I , J , K ) 

COF  4 ( I , J , K , 2 ) =AIM ( I , J, K) 

COF  4 ( I , J , K , 3 ) =AJP ( I , J , K ) 

COF  4 ( I , J , K , 4 ) =AJM ( I , J , K ) 

COF  4 ( I , J , K , 5 ) =AKP ( I , J, K ) 

COF  4 ( I , J , K , 6 ) =AKM ( I , J , K ) 

COF  4 ( I , J, K, 7 )=AP( I , J , K) 

703  CONTINUE 

IF( ITER.LE.l)  GO  TO  409 
DO  408  K=2 , N2 
DO  408  J=2 , M2 
DO  408  1=2, L2 

AP ( I , J , K ) =AP ( I , J, K ) /RELAX ( NP ) 

CON ( I , J , K ) =CON ( I , J , K )  +  ( 1 . -RE LAX ( NP  )  ) *AP( I,J,K)*P(I,J,K) 

408  CONTINUE 

409  CONTINUE 
CALL  TDMA 
NF=  1 
IST-3 
JST=2 
KST=2 

DO  704  1=1, LI 
DO  704  J=1 , Ml 
DO  704  K=1 , Nl 
AIP( I , J , K ) =COFl ( I , J, K  ,  1  ) 

AIM ( I , J , K )=COFl ( I , J , K , 2  ) 

AJP ( I , J,K)=COFl( I , J,K, 3) 


AJM( I , J , K ) =  COF 1 ( I , J , K, 4 ) 

AKP( I ,  J , K)=COFl( I ,J ,K, 5) 

AKM ( I  ,  J  , K )=C0F1 ( I , J , K, 6 ) 

AP(I,J,K)=C0Fl(I,J,K,7) 

CON ( I  ,  J  ,  K  ) =CONl ( I , J , K ) 

704  CONTINUE 

DO  413  K=2,N2 
DO  413  J=2 , M2 
DO  413  1=3 , L2 

413  CON ( I , J , K ) =CON ( I , J , K ) +DU ( I , J , K ) * AP ( I , J , K ) * ( P ( I - 1 , J , K ) — P ( I , J , K ) ) 
CALL  TDMA 
NF  =  2 
I  ST=2 
JST=3 
KST=2 

DO  714  1=1, Ll 
DO  714  J=1 , Ml 
DO  714  K=1,N1 
AIP( I ,J,K)=COF2( I ,J,K, 1 ) 

AIM( I , J , K ) =COF2 ( I , J , K , 2 ) 

AJP( I ,J,K)=COF2( I ,J,K, 3) 

AJM ( I , J, K)=COF2( I , J, K, 4 ) 

AKP(I,J,K)=COF2(I,J,K,5) 

AKM ( I , J , K ) = CO F2 ( I , J , K , 6 ) 

AP(I,J,K)=COF2(I,J,K,7) 

CON ( I , J , K ) =CON2 ( I , J , K ) 

714  CONTINUE 

DO  513  K=2,N2 
DO  513  1=2 , L2 
DO  513  J=3 , M2 

513  CON(I,J,K)=CON(I,J,K)+DV(I,J,K)*AP(I,J,K)MP(I,J-l,K)-P(I,J,K)  ) 
CALL  TDMA 
NF=3 
I  ST=2 
JST=2 
KST=  3 

DO  724  1=1, Ll 
DO  724  J=1 , Ml 
DO  724  K=l,Nl 
AIP(I,J,K)=COF3(I,J,K,l) 

AIM( I , J , K ) =COF 3 ( I , J, K,  2 ) 

AJP ( I , J,K)=COF3( I , J,K,  3) 

AJM ( I , J,K)=COF3( I , J , K  ,  4 ) 

AKP ( I , J , K ) -COF3 ( I , J , K , 5 ) 

AKM ( I , J , K ) =COF3 ( I , J , K ,  6 ) 

AP ( I,J,K)=COF3( I , J , K , 7 ) 

CON ( I, J,K)=CON3( I, J,K) 

724  CONTINUE 

DO  523  J=2,M2 
DO  523  1=2, L2 
DO  523  K=  3 , N2 

523  CON ( I , J , K ) =CON { I,J,K)+DW(I,J,K)*AP(I,J,K)*(P(I,J,K-1)-P(I,J,K)) 
CALL  TDMA 
C 

COEFFICIENTS  FOR  THE  PRESSURE  CORRECTION  EQUATION - 

C 

DO  706  1=1, Ll 
DO  706  J=1 , Ml 
DO  706  K=1 , Nl 
A I P ( I , J , K ) =COF 4 { I , J,  K,  1  ) 
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AIM( I  ,  J  ,  K  ) =C0F4 ( I , J, K, 2  ) 

AJP ( I ,  J  ,  K  ) =  C0F4 ( I , J, K, 3 ) 

AJM ( I,J,K)=C0F4( I , J , K , 4  ) 

AKP ( I , J,K) =C0F4 ( I , J, K, 5  ) 

AKM ( I , J, K ) =C0F4 ( I , J, K, 6 ) 

AP ( I  ,  J  ,  K  ) =C0F4 ( I , J, K, 7  ) 

706  CONTINUE 
NF=4 

I F ( . NOT . L SOLVE ( NF ) )  GO  TO  500 

I  ST=2 

JST=2 

KST=2 

CALL  DIFFUS 
SMAX=0 . 

SSUM=0  . 

DO  410  K=2 , N2 
DO  410  J=2,M2 
DO  410  1=2, L2 
VOL=YCV( J ) *XCV( I ) *ZCV( K ) 

410  CON (  I  ,  J  ,  K  ) =CON ( I , J , K ) * VOL 
DO  701  K=2,N2 
DO  701  J=2 , M2 

CON ( 2 , J , K ) =RHO ( 1 , J , K ) *U( 2 , J , K ) * YCV { J ) *ZCV( K ) 

701  CONTINUE 

DO  702  1=2 , L2 
DO  702  K=2,N2 

CON ( I , 2 , K)=RHO( I , 1 , K) *V( I , 2, K) *XCV( I ) *  ZCV ( K ) 

702  CONTINUE 

DO  153  1=2, L2 
DO  153  J=2 , M2 

CON ( I , J , 2 ) =RHO ( I,J,1)*W(I,J,2) *XCV( I ) *  YCV ( J ) 

153  CONTINUE 
C 

DO  351  K=2,N2 
DO  351  J=2 , M2 
DO  351  1=2, L2 
AREA= YCV ( J ) *ZCV( K ) 

ARHO-AREA* ( FX( 1  +  1 ) *RHO( 1  +  1 , J , K )  +  FXM ( 1  +  1  ) *RHO( I , J , K )  ) 
FLOW=ARHO*U( 1+1 , J,K) 

CON ( I, J,K) =CON ( I , J , K ) -FLOW 
CON ( 1+1, J,K)=CON( I + 1 , J , K ) +FLOW 
C 

AREA=XCV( I ) *ZCV ( K ) 

ARHO-AREAM  FY( J  +  l ) *RHO( T , J  +  l , K ) +FYM( J  +  l ) *RHO( I , J , K ) ) 
FLOW=ARHO*V( I , J+1,K) 

CON ( I , J , K ) =CON ( I , J , K ) -FLOW 
CON ( I , J+l , K ) =CON ( I , J+l , K ) +FLOW 
C 

AREA=XCV( I ) *YCV ( J ) 

ARHO-AREA* ( FZ ( K+l ) *RHO(  I , J , K  +  l ) +FZM( K  +  l ) *RHO(  I, J,K)  ) 
FLOW=ARHO*W( I , J , K+l ) 

CON ( I , J , K ) =CON ( I , J , K ) -FLOW 
CON ( I , J , K+ 1 ) =FLOW 
PC ( I , J , K ) =0 . 

351  CONTINUE 

DO  352  1=2, L2 
DO  352  J=2,M2 
DO  352  K  =  2 , N2 

SMAX=AMAXl ( SMAX , ABS ( CON ( I , J , K ) ) ) 

SSUM=SSUM+CON( I , J , K ) 
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352  CONTINUE 
CALL  TDMA 
C 

COME  HERE  TO  CORRECT  THE  VELOCITIES - 

C 

DO  511  K=2 , N2 
DO  511  J=2 , M2 
DO  511  1=2, L2 

I F ( I . NE . 2 )  U ( I , J , K ) =U { I , J , K ) +DU { I,J,K)*(PC(I-1,J,K)-PC(I,J,K)) 

I F ( J . NE . 2 )  V ( I , J , K ) =V ( I ,J,K)+DV( I , J , K ) * ( PC ( I  ,  J- 1  ,  K ) -PC  (  I  ,  J  ,  K )  ) 

I F ( K . NE . 2 )  W(I,J,K)=W(I,J,K) +  DW (I,J,K)*(PC(I,J,K-1)-PC(I/J,K) ) 
511  CONTINUE 
500  CONTINUE 
C 

COEFFICIENTS  FOR  OTHER  EQUATIONS - 

C 

CALL  RESET 
I  ST=2 
J  ST=  2 
KST=2 

DO  600  NF=5,NFMAX 

IF( . NOT . LSOLVE ( NF ) )  GO  TO  600 

CALL  DIFFUS 

REL= 1 . -RELAX ( NF ) 

DO  601  1=2, L2 
DO  601  J=2 , M2 
DO  601  K=2 , N2 

COEFFICIENTS  AWEST  AND  AEAST 
AREA=  YCV ( J ) *  ZCV ( K ) 

FLOW=AREA  *  U ( I+1,J,K)*(FX(I+1) *RHO( I + 1 , J , K ) +FXM (  1  +  1 ) *RHO(  I , J , K )  ) 
DIFF=AREA*2 . *GAM ( I , J , K ) *GAM ( 1  +  1 , J , K )/( XCV ( I ) *GAM ( 1+1 , J, K)+ 

1  XCV (  I  +  1 ) A  GAM (I,J,K)+1.0E-30) 

CALL  PROFIL 

AI M ( I + 1 , J , K ) =ACOF+AMAXl ( ZERO, FLOW) 

AIP( I, J,K)=AIM( 1+1, J,K)-FLOW 
COEFFICIENTS  ANORTH  AND  ASOUTH 
AREA=  XCV ( I ) *  ZCV ( K ) 

FLOW=AREA* V ( I , J+ 1 , K ) * ( FY( J+l ) *RHO( I , J  +  l , K ) +  FYM ( J+ 1 ) *RHO( I , J, K )  ) 
D I F  F =AREA*  2 . *GAM ( I , J , K ) *GAM ( I , J+l , K )/( YCV ( J ) *GAM ( I , J  +  l , K  )  + 

1  YCV( J+l ) *GAM( I , J , K ) +1 . 0E-30 ) 

CALL  PROFIL 

AJM ( I , J+l,K)=ACOF+AMAXl( ZERO, FLOW) 

AJ  P ( I, J,K)=AJM( I , J  + 1 , K ) - FLOW 
COEFFICIENTS  AOUT  AND  AINTO 
AREA=  YCV( J ) *XCV( I ) 

FLOW=AREA*W ( I , J , K+ 1 ) * ( FZ ( K+ 1 ) *RHO( I , J, K+l )+FZM( K+l ) *RHO( I , J, K) ) 
DI FF-AREA*  2 . *GAM ( I , J, K) *GAM( I , J, K  +  l )/( ZCV ( K ) *GAM ( I , J, K+l )  + 

1  ZCV ( K+l ) *GAM ( I , J , K ) + 1 . 0E-30 ) 

CALL  PROFIL 

AKM ( I , J , K  +  l ) =ACOF +AMAX1 ( ZERO, FLOW) 

AKP ( I , J , K ) =AKM ( I , J , K+l ) -FLOW 
601  CONTINUE 

DO  610  J=2,M2 
DO  610  K=2 , N2 

COEFFICIENTS  AWEST  AND  AEAST 
AREA=YCV ( J ) *ZCV( K ) 

FLOW=AREA*U ( 2 , J , K ) *RHO( 1 , J , K) 

D I F  F =AREA  *GAM ( 1 , J , K )/XDIF( 2 ) 

C  A r  •  L  PROFIL 

A I M ( 2 , J , K ) =ACOF  +  AMAXl ( ZERO, FLOW) 
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FLOW=AREA*U  (  Ll  ,  J  ,  K)  *RHO(  ’-1  ,  J,  K  ) 

DIFF=AREA*GAM( Ll , J,K)/XDIF(Ll ) 

CALL  PROFIL 

AIP( L2 , J , K ) =ACOF +AMAX1 ( ZERO , FLOW) -FLOW 

610  CONTINUE 

DO  611  1=2, L2 
DO  611  K=2,N2 

COEFFICIENTS  ANORTH  AND  ASOUTH 
AREA=  XCV ( I ) *  ZCV ( K ) 

FLOW=AREA* V ( I , 2 , K ) *RHO( I , 1 , K ) 

DIFF=AREA*GAM( 1,1, K)/YDIF( 2 ) 

CALL  PROFIL 

A JM ( I , 2 , K ) =ACOF+AMAXl ( ZERO, FLOW) 

FLOW=AREA* V ( I , Ml , K ) *  RHO ( I , Ml , K ) 

DIFF=AREA*GAM( I , Ml , K ) /YDI F ( Ml ) 

CALL  PROFIL 

AJ  P ( I , M2 , K )=ACOF+AMAXl ( ZERO, FLOW) -FLOW 

611  CONTINUE 

DO  612  1=2, L2 
DO  612  J=2 , M2 

COEFFICIENTS  AOUT  AND  AINTO 
AREA=  YCV ( J ) *XCV ( I ) 

FLOW=AREA*W( I , J , 2 ) *RHO( I , J , 1 ) 

DIFF=AREA*GAM( I,J,1)/ZDIF(2) 

CALL  PROFIL 

AKM ( I , J , 2 ) =ACOF+AMAXl ( ZERO, FLOW) 

FLOW=AREA* W ( I , J , Nl ) *RHO( I , J , N1 ) 

DIFF=AREA*GAM( I , J , Nl ) /ZDI F ( Nl ) 

CALL  PROFIL 

AKP ( I , J,N2 )=ACOF+AMAXl( ZERO, FLOW) -FLOW 

612  CONTINUE 

DO  3987  1=1, Ll 
DO  3987  J=1 , Ml 
DO  3987  K= 1 , Nl 
VOL  =  YCV ( J ) *XCV( I ) *  ZCV ( K ) 

APT=RHO ( I , J , K)/DT 

AP ( I , J , K ) =AP ( I , J , K ) -APT 

CON ( I , J , K)=CON( I , J , K ) +APT*  F ( I , J , K,NF ) 

AP ( I , J , K ) 

1  =( -AP( I , J , K ) * VOL+AI P ( I , J , K ) +AI M ( I , J , K)+AJP( I , J , K ) +AJM ( I , J, K) 

2  +AKM( I ,J,K)+AKP( I ,J,K) ) 

3/RELAX ( NF ) 

CON ( I , J , K ) =CON ( I , J,K) *VOL+REL*AP( I , J, K) *F( I , J , K , NF ) 

C  PRINT* , I , J , K 

C  PRINT* , AIM( I , J , K ) , A JM ( I , J , K ) ,AKM(I,J,K) 

C  PRINT* ,AIP( I , J, K) , AJP( I , J, K) ,AKP( I , J , K ) 

C  PRINT* , AP( I , J , K )  , CON( I , J , K  ) 

3987  CONTINUE 
CALL  TDMA 
600  CONTINUE 

TI ME=TI ME+DT 
I TER= I TER+ 1 

I F ( ITER. GE. LAST)  LSTOP=.TRUE. 

RETURN 

END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SUBROUTINE  OPTION 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
LOGICAL  LSOLVE, LPRINT, LBLK, LSTOP 

COMMON  F(17,17,13,5),P(17,17,13)  , RHO ( 1 7 , 1 7 , 1 3 ) , GAM ( 1 7 , 1 7 , 1 3  )  , 


M 


o  o  n  o 


1  CON ( 17,17,13) , AKP( 17,17,13) , AKM( 17 , 17, 13), AP (17, 17, 13), 

2  AI P ( 1 7 , 1 7 , 1 3 ) , AIM ( 17 , 1 7 , 1 3 ) , AJP ( 17 , 1 7 , 1 3 ) ,AJM( 17,17,13) 

COMMON 

3  X ( 1 7 )  , XU (17),XDIF(17) ,XCV( 17 ) ,XCVS( 17  )  , 

4  Y( 17) , YV( 17 ) , YDIF( 17) , YCV( 17 ) , YCVS( 17 ) , 

5  Z( 17 ) ,ZW( 17 ) ,ZDIF( 17 ) , ZCV(17) , ZCVS(17) , 

6  YCVR ( 17 ) , YCVRS( 17 ) , ARX( 17 ) , ARXJ ( 17 ) , ARXJP( 17 ) , 

7  R ( 1 7 ) , RMN (17) ,SX(17) , SXMN(17) ,XCVI(17) ,XCVIP(17) , 

8  YCVJ ( 17 ) , YCVJP ( 17 ) , ZCVK ( 17 ) , ZCVKP( 17 ) 

COMMON  DU (17, 17, 13), DV (17, 17, 13), DW (17, 17, 13), FV (17), FVP( 17 ) , 

1  FX( 17 )  , FXM ( 17 ) , FY( 17 ) , FYM ( 1 7 ) , PT ( 1 7 ) , QT ( 1 7 )  , 

2  FZ( 17 ) , FZM( 17  ) 

COMMON/I NDX/NF , NFMAX , NP , NRHO , NGAM ,Ll,L2,L3,Ml,M2,M3,Nl,N2,N3, 

1 1  ST , JST , KST , ITER, LAST, TITLE ( 13 )  , RELAX ( 1 3 )  , TI ME , DT , XL , YL , Z L , 
2IPREF, JPREF, KPREF, LSOLVE( 10 )  , LPR I NT ( 13 )  ,  LBLK(  11 )  , MODE , NTIMES ( 10 )  , 
3RHOCON, ZERO 

COMMON/ 1 NDX/RELAX( 13)  ,LPRINT(13)  ,LBLK(11)  , NTIMES  (  10 )  , 

1 LSOLVE ( 1 0 ) , TIME , DT ,XL,YL,ZL,S, RHOCON, ZERO, 

2NF, NFMAX, NP, NRHO, NGAM, Ll,L2,L3,Ml,M2,M3,Nl,N2,N3, 

3  1ST, JST, KST, ITER, LAST, 

4  IP REF, JPREF, KPREF, MODE 
COMMON/HEADIN/TITLE 
CHARACTER* 1 0  TITLE(13) 

COMMON/CNTL/LSTOP 
COMMON/SORC/SMAX , SSUM 
COMMON/COEF/FLOW, DIFF, ACOF 

DIMENSION  U(17,17,13),V(17,17,13),W(17,17,13),PC(17,17,13) 
DIMENSION  T  (  1 7 , 1 7 , 1 3  ) 

EQUIVALENCE ( F ( 1 , 1 , 1 , 1 ) ,U(1,1,1)),(F(1,1,1,2),V(1,1,1)), 

1  (F(l,  1,1,3  )  ,W( 1,1,1)  )  ,(F(1,1,1,4)  , PC (1,1,1)  ) 

2  , (F(l,l,l,5) ,T(1,1,1)  ) 


10  FORMAT ( 26(1H*),3X,A10,3X,26(1H*)) 

20  FORMAT(lX,4H  I  =, 16,619) 

30  FORMAT( IX, 1HJ ) 

40  FORMAT( IX, 12 , 3X, IP, 7E9 . 2  ) 

50  FORMAT ( 1 H  ) 

51  FORMAT ( IX , ' I  =' , 2X, 7 ( 14 , 5X) ) 

52  FORMAT ( IX , ' X  =',1P,7E9.2) 

53  FORMAT ( ' TH  =',1P,7E9.2) 

54  FORMAT (IX, ' J  = '  , 2X , 7 ( I  4 , 5X )  ) 

55  FORMAT ( IX , ' Y  =',lP,7E9.2) 

56  FORMAT ( IX , ' K  = ' , 2X , 7 { I  4 , 5X ) ) 

57  FORMAT(lX,'Z  =',1P,7E9.2) 

59  FORMAT (IX, ' K  =',2X,I4) 

C*****************************************  ****** 

ENTRY  UMESH 
XU( 2 ) =0 . 

DX=XL/FLOAT( Ll-2  ) 

DO  1  1=3, LI 

1  XU ( I ) =XU( 1-1 ) +DX 
YV ( 2  )  =0  . 

DY=YL/ FLOAT ( Ml - 2 ) 

DO  2  J=  3 , Ml 

2  YV( J )=YV( J-l ) +DY 
ZW( 2 ) =0  . 

DZ=ZL/FLOAT( Nl-2 ) 

DO  3  K=  3 , N1 

3  ZW ( K ) =ZW ( K  —  1 ) +DZ 


RETURN 


ENTRY  PRINT 


C**********************] 


IF(  .NOT.LPRINT{ 3 ) )  GO  TO  80 
80  CONTINUE 


C 


I F ( .NOT. LPRINT( NP ) )  GO  TO  90 
C 

CONSTRUCT  BOUNDARY  PRESSURES  BY  EXTRAPOLATION 
C 

DO  91  K=2 , N2 
DO  91  J=2 , M2 

P( 1,J,K)=(P(2,J,K) *XCVS ( 3)-P( 3fJ,K)*XDIF(2) )/XDIF( 3 ) 

91  P(L1,J,K)=(P(L2,J,K) *XCVS (L2)-P(L3,J,K)*XDIF(L1) )/XDIF(L2) 
DO  92  K=2 , N2 

DO  92  1=2, L2 

P(I,1,K)=(P(I,2,K) * YCVS ( 3)-P(I,3,K)*YDIF(2) )/YDIF( 3 ) 

92  P( I , Ml , K ) = ( P ( I ,M2 , K ) * YCVS { M2 ) -P ( I , M3 , K ) * YDI F ( Ml ) )/YDIF(M2 ) 
DO  93  J=2 , m2 

DO  93  1=2, L2 

P(I,J,1)=(P(I,J,2) *ZCVS( 3)-P(I,J,3)*ZDIF(2) )/ZDIF(3) 

93  P(I,J,Nl)  =  (P(I,J,N2) *  ZCVS ( N2 )-P(I,J,N3)*ZDIF(Nl))/ZDIF(N2) 
DO  94  K=2,N2 

P(1,1,K)=P(2,1,K)+P(1,2,K)-P(2,2,K) 

94  CONTINUE 

DO  95  J=2 , m2 

P(1,J,1)=P(2,J,1)+P(1,J,2)-P(2,J,2) 

95  CONTINUE 

DO  96  1=2, L2 

P(I,1,1)=P(I,2,1)+P(I,1,2)-P(I,2,2) 

96  CONTINUE 

P(l,l,l)=(P(l,l,2)+P(l,2,l)+P(2,l,l))/3.0 
P( LI , 1 , 1 )  =  ( P( L2,l , 1 )+P<  LI ,2, 1 )+P(Ll , 1, 2) )/3. 0 
P(l,l,Nl)=(P(l,l,N2)+P(l,2,Nl)+P(2,l,Nl))/3.0 
P(1,M1,1)=(P(1,M2,1)+P(1,M1,2)+P(2,M1,1) )/3.0 
P( LI ,M1 , 1 )=( P( L2,M1 , 1 )+P{ LI ,M2, 1 )+P< LI ,Ml , 2 ) )/3 . 0 
P( 1 ,Ml , N1 )=( P( 2 ,Ml ,Nl ) +P( 1 , M2 ,N1 )+P( 1 ,M1 ,N2 ) )/3 . 0 
P( LI ,M1 ,N1 )=(P( L2,Ml ,Nl )+P{ LI ,M2 ,N1 )+P( LI ,M1 ,N2 ) )/3 . 0 
PREF  =P ( IPREF, JPREF, KPREF) 

DO  97  K= 1 , Nl 
DO  97  J=1 , Ml 
DO  97  1=1, LI 

97  P( I , J, K)=P( I , J, K )-PREF 
90  CONTINUE 

C 

PRINT  50 
I END=0 

301  I F ( IEND.EQ.L1 )  GO  TO  310 
IBEG=IEND+1 

I END= I END  +  7 
I END=MI NO ( I  END , Ll  ) 

PRINT  50 

PRINT  51 , ( I , I=IBEG, IEND) 

IF( MODE. EQ. 3 )  GO  TO  302 
PRINT  52, (X( I ), I-IBEG, IEND) 

GO  TO  303 

302  PRINT  53, (X( I ), I-IBEG, IEND) 

303  GO  TO  301 
310  J END=  0 


Yj 


PRINT  50 

311  I F { JEND.EQ.M1 )  GO  TO  320 
JBEG=JEND+1 
J  END=JEND+7 
JEND=MINO(JEND,Ml  ) 

PRINT  50 

PRINT  54,  (  J, J  =  JBEG, JEND) 

PRINT  55, (Y( J) , J=JBEG, JEND) 

GO  TO  311 

320  KEND=0 
PRINT  50 

321  IF( KEND. EQ.Nl )  GO  TO  330 
KBEG=KEND+ 1 
KEND=KEND+7 
KEND=MIN0 ( KEND , Nl ) 

PRINT  50 

PRINT  56, (K,K=KBEG, KEND ) 

PRINT  57, (Z(K) , K=KBEG, KEND) 

GO  TO  321 
330  CONTINUE 
C 

DO  999  NF=1 , NGAM 
I F  (  . NOT . LPRINT ( NF ) )  GO  TO  999 
PRINT  50 

PRINT  10 , TITLE( NF ) 

DO  998  K=1 , Nl 
PRINT  50 
PRINT  59, K 
I FST= 1 
J  FST=1 
KFST=1 

I F  (  NF . EQ . 1 . OR . NF . EQ . 4 )  IFST=2 
IF(NF.EQ.2.0R.NF.EQ.4)  JFST=2 
IF(NF.EQ.3.0R.NF.EQ.4)  KFST=2 
IBEG=I FST-7 
110  CONTINUE 

IBEG=IBEG+7 
IEND=IBEG+6 
I END=MIN0 ( I  END , Ll ) 

PRINT  50 

PRINT  20 , ( I , I=IBEG, IEND ) 

PRINT  30 
JFL-JFST+Ml 
DO  115  J J=JFST , Ml 
J=JFL-J J 

PRINT  40, J,  (F(I,J,K,NF)  , I  =  I  BEG , I  END ) 

115  CONTINUE 

I F ( IEND.LT.Ll)  GO  TO  110 

998  CONTINUE 

999  CONTINUE 
RETURN 
END 

C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

C  PROBLEM  DEPENDENT  PORTION  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

BLOCK  DATA 

LOGICAL  LSOLVE , LPRINT, LBLK, LSTOP 

COMMON  F(17,17,13,5),P(17,17,13) , RHO ( 1 7 , 1 7 , 1 3 ) , GAM ( 1 7 , 1 7 , 1 3 ) , 

1  CON ( 17,17,13) ,AKP( 17,17, 13) ,AKM( 17, 17,13) ,AP(17,17,13), 
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2  AIP( 17,17, 13) ,AIM( 17,17,13) ,AJP( 17 ,17, 13) ,AJM( 17,17,13) 

COMMON 

3  X ( 1 7 ) , XU (17),XDIF(17) ,XCV( 17 ) ,XCVS( 17 ) , 

4  Y( 17 ) , YV( 17 ) , YDIF( 17 ) , YCV( 17 ) , YCVS( 17 ) , 

5  Z ( 1 7 ) , ZW{ 1 7 ) , ZDI F ( 17 ) ,ZCV(17) ,ZCVS(17) , 

6  YCVR ( 17 ) , YCVRS (17) , ARX( 17 ) , ARXJ ( 17 ) , ARXJP( 17 ) , 

7  R ( 1 7 )  , RMN (17) , SX (  1 7 )  ,SXMN(17) ,XCVI(17) ,XCVIP(17) , 

8  YCVJ( 17 )  , YCVJ P ( 17 ) , ZCVK( 17 ) , ZCVKP( 17 ) 

COMMON  DU (17, 17, 13), DV (17, 17, 13), DW (17, 17, 13), FV (17), FVP( 17 ) , 

1  FX( 17 ) , FXM ( 17 ) ,FY( 17 ) , F YM ( 17 ) , PT( 17 ) , QT( 1 7 ) , 

2  FZ( 17 ) , FZM ( 17  ) 

COMMON/ I NDX/RE LAX ( 13), LPRINT{ 13), LBLK ( 11 ) ,NTIMES( 10 ) , 

1 LSOLVE ( 10) , TI ME , DT ,XL,YL,ZL, S, RHOCON, ZERO, 

2NF , NFMAX , NP , NRHO , NGAM ,Ll,L2,L3,Ml,M2,M3,Nl,N2,N3, 

3 1  ST, JST, KST, ITER, LAST, 

4IPREF, JPREF, KPREF, MODE 
COMMON/HEAD IN/TITLE 
CHARACTER* 10  TITLE(13) 

COMMON/CNTL/LSTOP 
COMMON/SORC/SMAX , SSUM 
COMMON/COEF/FLOW, DIFF, ACOF 

DIMENSION  U(  1 7 , 1 7 , 1 3  ) ,V(17,17,13) ,W( 17,17,13) ,PC(17,17,13) 
DIMENSION  T( 17 , 17 , 13 ) 

EQUIVALENCE(F( 1,1, 1,1)  ,U(  1,1,  1)  ),(F(1,1,1,2),V(1,1,1)), 

1  (F( 1,1, 1,3) ,W( 1,1,1)  )  ,  (F( 1,1, 1,4)  , PC (1,1,1)  ) 

2  , ( F( 1 , 1 , 1 , 5 ) ,T( 1 , 1 , 1 )  ) 

DATA  NFMAX , NP , NRHO , NGAM/5 ,6,7, 8/ 

DATA  LSTOP, LSOLVE, LPRINT/24* . FALSE./ 

DATA  MODE, TIME, ITER/1,0. ,0/ 

DATA  RELAX, NTIMES/13*1 . 0 , 10*1/ 

DATA  LBLK/11* .TRUE./ 

DATA  DT, IPREF, JPREF, KPREF, RHOCON/1 . E+ 1 0 , 1  ,  1  ,  1  ,  1 . 0/ 

C  NF=1 , 2 , 3  STAND  FOR  U,V  AND  W  VELOCITIES. 

C  NF-6  IS  FOR  PRESSURE 
C  NF=  5  IS  FOR  TEMPERATURE 

C  L SOL VE= TRUE  SOLVES  THAT  PARTICULAR  PHI 
DATA  ( LSOLVE( I ) , 1=1 , 6 )/6* . TRUE./ 

C  LPRINT( NF ) =TRUE  PRINTS  VARIABLE  ASSOCIATED  WITH  NF  ON  CALLING  PRINT 
DATA  LPRINT( 5)/l* .TRUE./ 

C  TERMINATE  ITERATIONS  AT  I TER=  LAST 
DATA  LAST/100/ 

C  UNDERELAXATION  FACTORS 

DATA  RELAX ( 1 ) , RELAX ( 2 ) , RELAX ( 3 ) , RELAX ( 4) /0. 5, 0.5, 0.5, 1.0/ 

DATA  RELAX( 5 ) , RELAX( 6 )/l .0,1. 0/ 

C  TITLES  FOR  THE  FIELD  PRINTOUTS 

DATA  TITLE ( 1 ) , TITLE ( 2 )  , TITLE ( 3 ) , TITLE ( 5 ) , TITLE  (  6 ) / 
l'U' , 'V' , 'W' , 'T' , 'P '/ 

C  NUMBER  OF  SWEEPS  IN  THE  LINE-BY-LINE  TDMA  ALGORITHM 
DATA  NTI MES ( 3) ,NTIMES( 4 ) /2  *  5/ 

DATA  ZERO/0./ 

END 

(]************************************************************ 
SUBROUTINE  USE 

C******A***************************************************** 

C  IF  LARGER  NUMBER  OF  GRID  POINTS  IS  TO  BE  USED,  THE  DIMENSION 
C  STATEMENTS  MUST  BE  CHANGED  THROUGH  THE  PROGRAM  TO  ACCOMODATE 
C  VALUES  GRAETER  THAN  (17,17,13)  ETC. 

C************************************************************ 

LOGICAL  LREAD , LWRITE 


4  L 


n  o  no 


LOGICAL  LSOLVE, LPRINT, LBLK, LSTOP 

COMMON  F(17,17,13,5),P(17,17,13) ,RHO( 17 , 17,13) ,GAM( 17, 17,13), 

1  CON ( 17,17,13) ,AKP(17,17,13) , AKM ( 17, 17, 13), AP (17, 17, 13), 

2  AIP( 17,17,13)  ,AIM( 17,17,13) , A J P ( 1 7 , 1 7 , 1 3 ) ,AJM(17,17,13) 

COMMON 

3  X ( 1 7 ) , XU (17),XDIF(17) ,XCV( 17 ) ,XCVS{ 17 ) , 

4  Y ( 1 7 )  , YV ( 17 ) , YDI F ( 1 7 ) ,YCV(17) ,YCVS(17)  , 

5  Z(17),ZW(17),ZDIF(17) ,ZCV(17) ,ZCVS(17)  , 

6  YCVR( 17 ) , YCVRS ( 17 ) , ARX( 17 ) , ARXJ( 17 ) , ARXJP( 17 ) , 

7  R( 17 ) , RMN ( 17 ) , SX( 17 ) , SXMN ( 17 ) ,XCVI ( 17 ) ,XCVIP( 17 ) , 

8  YCVJ ( 17 ) , YCVJP ( 17 ) , ZCVK( 17 ) , ZCVKP ( 17 ) 

COMMON  DU( 17,17,13) ,DV(17,17,13) , DW (17,17,13),FV(17), FVP(  17 )  , 

1  FX( 17 ) , FXM( 17 ) , FY( 17 )  , FYM(  17 )  , PT ( 1 7 ) , QT ( 1 7 ) , 

2  FZ ( 1 7 ) , FZM( 17 ) 

COMMON/ I NDX/RELAX ( 13 ) , LPRINT ( 13 ) , LBLK ( 11 ) , NTIMES( 10 ) , 

1 LSOLVE ( 10) ,TIME,DT,XL,YL,ZL,S, RHOCON, ZERO, 

2NF , NFMAX , NP , NRHO , NGAM ,Ll,L2,L3,Ml,M2,M3,Nl,N2,N3, 

3  I  ST, JST, KST, ITER, LAST, 

4  I PREF, JPREF, KPREF , MODE 
COMMON/HEAD IN/TITLE 
CHARACTER* 1 0  TITLE(13) 

COMMON/CNTL/ LSTOP 
COMMON/SORC/SMAX , SSUM 
COMMON/COEF/FLOW , DI FF , ACOF 

DIMENSION  U ( 1 7 , 1 7 , 1 3 ) ,V(17,17,13) ,W(17,17,13) ,PC(17,17,13) 
DIMENSION  T( 17, 17, 13) , BUOY ( 1 7 , 1 7 , 1 3  ) 

C  SPECIFY  EQUIVALENCE  FOR  QUANTITIES  INSIDE  AND  OUTSIDE  SUBROUTINE  USE 
C  T=TEMPERATURE, PC=PRESSURE  CORRECTION 

EQUIVALENCE ( F ( 1 ,1,1,1),U(1,1,1)),(F(1,1,1,2),V(1,1,1)), 

1  (F(1,1,1,3),W(1,1,1)),(F(1,1,1,4), PC (1,1,1)) 

2  , ( F( 1 , 1 , 1 , 5 ) ,T( 1 , 1 , 1 )  ) 

DIMENSION  ANUZ( 13) ,ANU3PT( 13) 

CHARACTER* 1  AREAD 

*********************************************** 

ENTRY  MESH 

DOMAIN  LENGHTS  IN  THE  3  DIRECTIONS 
XL=1 . 0 
YL= 1  .  0 
ZL=1 . 0 

C  NUMBER  OF  GRID  POINTS  IN  THE  3  DIRECTIONS 
Ll  =  17 
Ml-17 
Nl  =  l  3 

C  INVOKE  THE  UNIFORM  MESH  OPTION 
CALL  UMESH 
RETURN 
C 

ENTRY  BEGIN 

C  ITERATIONS  STOP  AT  I TER=LAST 

WRITE( * , * ) ' NO.  ITERATIONS' 

READ ( * , * ) LAST 
C  READ  RAYLEIGH  NUMBER 
WR I TE ( * , * ) ' RALI= ' 

READ ( * , * ) RA 
C  READ  PRNDTL  NUMBER 

WR I TE ( * , * ) ' PRANTL= ' 

READ ( * , * ) PR 

C  DETERMINE  IF  WANT  TO  USE  A  PREVIOUSLY  COMPUTED  SOLUTION  AS 
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n  o 


C  AN  INITIAL  GUESS 

WRITE ( * , * )  ' READ  FROM  INPUT  FILE  (0/1)' 

READ (  * ,  *  ) I  READ 
I F (  I  READ . EQ . 1 ) LREAD= . TRUE . 

I F ( I  READ. EQ. 0 ) LREAD= . FALSE . 

C  PROVIDE  INITIAL  GUESS.  THE  PROGRAM  SOLVES  FOR  THE  INTERIOR  POINTS 
C  ONLY.  HENCE  THE  BOUNDARY  VALUSE  FOR  THE  TEMPERATURES  AT  THE  HOT 
C  AND  COLD  BOUNDARIES  HAVE  BEEN  ALREADY  SPECIFIED. 

DO  101  1=1, LI 
DO  101  J=1 , Ml 
DO  101  K=1 , Nl 
U ( I , J , K ) =0  . 

V  (  I  ,  J  ,  K  )  =  0  . 

W  (  I  ,  J , K ) =0 . 

T( I , J,K)»  0. 

C  HOT  WALL  AT  X=0. 

T(1,J,K)=1.0 
101  CONTINUE 

I F ( . NOT . LREAD )  RETURN 
C  READ  DATA  FROM  INPUT  FILE 
REWIND  (4) 

READ ( 4 )  X, Y,Z,XU,YV,ZW,U,V,W, P,T 
RETURN 

★  ★A********'************************************ 

ENTRY  VARRHO 

(]*****HH*****H*******H******H*********H*** 

RETURN 

C 

Q**H*H****************H***HH*************** 

C  INCORPORATE  BOUNDARY  CONDITIONS 
ENTRY  BNDRY 
DO  864  1=1, LI 
DO  864  J= 1 , Ml 

C  UPDATE  VELOCITIES  AND  TEMPERATURE  AT  THE  SYMMETRY  LINE  (Z=ZL) 

C  ACORDING  TO  THE  ZERO  GRADIENT  BOUNDARY  CONDITION  FOR  U,V  AND  T 
U( I, J,Nl )=U( I , J,N2) 

V(  I  ,  J , Nl ) = V ( I , J , N2  ) 

T(  I , J , Nl )=T( I , J , N2  ) 

C  ADIABATIC  SIDE  (Z=0) 

T ( I , J , 1 ) =T ( I , J , 2 ) 

864  CONTINUE 

DO  865  1=1, LI 
DO  865  K=1 , Nl 
C  ADIABATIC  BOTTOM  (Y=0) 

T  (  I  ,  1 ,  K ) =T ( 1 , 2 , K ) 

C  ADIABATIC  TOP  (Y=YL) 

T( I,Ml,K)=T( I , M2 , K ) 

865  CONTINUE 

RETURN 

C 

ENTRY  PRTOUT 
IF( ITER.NE.0)  GO  TO  400 
PRINT  401 

401  FORMAT  ( IX, 'SIMPLER' , // ) 

400  CONTINUE 

C  MONITOR  SSUM , SMAX  AND  OTHER  QUANTITES  AS  ITERATIONS  PROCEED 
C  ON  CONVERGENCE,  SMAX  SHOULD  BE  VERY  SMALL  (LESS  THAN  1.0E-04) 

C  SSUM  SHOLD  ACHIEVE  A  SMALL  VALUE  WITHIN  A  FEW  ITERATIONS,  WELL 

C  BEFORE  CONVERGENCE.  SSUM  WILL  NOT  BE  SMALL  IF  THE  BOUNDARY 


C  CONDITIONS  ARE  NOT  WRITTEN  CORRECTLY. 

PRINT  403 , ITER, SSUM, SMAX, U( 5,5,9),V(5,5,5),T(8,6,5) 

403  FORMAT  (I6,5E12.4) 

COMPUTE  LOCAL  AND  GLOBAL  NUSSELT  NUMBERS  AS  ITERATIONS  PROCEED 
I F (  ( ITER . EQ. 250 ) 

2 .OR. ( ITER. EQ. 500  ) 

3. OR. ( ITER. EQ. 750) 

4 .OR. ( ITER. EQ. 1000  ) 

5 .OR. ( ITER. EQ. 1250 ) 

6 .OR. ( ITER. EQ. 1500  ) 

7 .OR. ( ITER. EQ. 1750 ) 

8. OR. ( ITER. EQ. 2000) 

1 .OR. ( ITER. EQ. 2250  ) 

2 .OR. ( ITER. EQ. 2500  ) 

3. OR. ( ITER. EQ. 2750) 

4 .OR.  ( ITER. EQ. 3000  ) 

5 .OR. ( ITER. EQ. LAST) ) THEN 

WRITE  ( 9 , * ) ' ITER, SSUM, SMAX, U( 5,5,9),V(5,5,5),T(8,6,5)' 

WRITE  ( 9 , * ) ITER, SSUM, SMAX, U( 5,5,9),V(5,5,5),T(8,6,5) 

COMPUTE  AVERAGE  NUSSELT  NUMBER 
ANUHOT=0 . 0 
ANUCLD=0 . 0 
DO  666  K=1,N1 
ANUZ ( K ) =  0 . 0 
DO  667  J  =  1 , M 1 

COMPUTE  AVERAGE  NUSSELT  NUMBER  AT  THE  HOT  AND  COLD  WALL 

ANUHOT=ANUHOT+ ( T ( 1,J,K)-T(2,J,K)) *YCV( J ) * ZCV ( K ) / ( XDI F ( 2 ) *YL*ZL ) 
ANUCLD=ANUCLD+(T(M2, J,K)-T(Ml , J,K)  ) *  YCV { J ) * ZCV ( K ) / ( XDI F ( Ll ) *YL*ZL) 
COMPUTE  LOCAL  NUSSELT  NUMBER  AT  THE  HOT  WALL 

ANUZ(K)=ANUZ(K)  +  (T(1, J,K)-T(2, J,K)  ) * YCV ( J ) / ( XDI F ( 2 ) * YL ) 

1  *  YCV (J)/(XDIF(2)*YL) 

667  CONTINUE 

C  WRITE  NUSSELT  NUMBERS  TO  AN  OUTPUT  FILE 

WRITE(9,*)'Z(K)=',Z(K), ' ANUZ ( K ) = ' , ANUZ ( K ) 

666  CONTINUE 

WRITE ( 9 , * ) ' ANUHOT  =', ANUHOT 
WRITE( 9 , * ) ' ANUCLD  =',ANUCLD 

I F ( ITER. EQ. LAST)  WRITE ( 9 , * ) ' RA= ' , RA , ' PR= ' , PR , ’ Z L- ' , Z L 
WRITE (9, *)'***************************************' 

ENDIF 

IF( ITER. NE. LAST)  RETURN 

C  GET  A  FIELD  PRINTOUT  AFTER  ITERATIONS  STOP 
CALL  PRINT 

C  WRITE  IN  ORDER  TO  COMMENCE  NEW  PROBLEM 
REWIND ( 4 ) 

WRITE ( 4)  X , Y , Z , XU , YV , ZW , U , V , W , P , T 
RETURN 

£AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

ENTRY  DIFFUS 
IF  (NF.EQ.4)  RETURN 
DO  500  1=1, Ll 
DO  500  J=1 , Ml 
DO  500  K-l.Nl 

C  DIFFUSIVITY  FOR  THE  U,  V  OR  W  EQUATIONS 
GAM ( I , J, K )=PR 

C  SPECIFY  ZERO  DIFFUSIVITY  FOR  THE  ZERO  GRADIENT  CONDITION  AT  Z=ZL 
C  FOR  THE  U  AND  V  VELOCITIES. 

I F ( NF.NE. 3 )GAM( I , J , Nl ) =0 . 0 
C  DIFFUSIVITY  FOR  THE  ENERGY  EQUATION 
IF  (NF.EQ.5)  THEN 
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GAM ( I, J,K)»1.0 

C  SPECIFY  ZERO  DIFFUSIVITIES  FOR  THE  ADIABATIC  BOUNDARIES 
GAM ( I  ,  1  ,  K ) =0 . 0 
GAM( I , Ml , K ) =0 . 0 
GAM( I , J , Nl ) =0 . 0 
GAM ( I , J, 1 )=0 . 0 
END  I  F 

500  CONTINUE 

IF  (NF.NE.2)  RETURN 

C  SOURCE  TERMS  ARE  EVALUATED  ONLY  FOR  INTERIOR  POINTS 
DO  501  1=2, L2 
DO  501  J=  3 , M2 
DO  501  K=2 , N2 

C  INTERPOLATE  TO  GET  THE  VALUE  OF  TEMPERATURE  AT  THE  CONTROL  VOLUME 
C  INTERFACE,  SINCE  THE  VELOCITY  U  IS  EVALUATED  AT  THE  INTERFACE. 

TM=  FY ( J ) *T( I , J,K)+FYM( J)*T( I, J-1,K) 

CONl=RA*PR*TM 

C  SUPPLY  ONLY  PART  OF  THE  SOURCE  TERM  IN  THE  MOMENTUM  EQ .  TO  AVOID 
C  DIVERGENCE  BEFORE  200  ITERATIONS 

CON ( I , J , K )= FLOAT ( ITER*  *2 ) *CONl/2  00 . **2 
IF  ( ITER.GT.200)CON( I, J,K)=CONl 

501  CONTINUE 

RETURN 

END 
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APPENDIX  C  -  Natural  convection  in  a  rectangular  box. 
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The  problem  being  considered  is  that  of  natural  convection  in  a  fluid- 
filled  rectangular  box.  Two  vertical  opposite  walls  are  maintained  at  a  hot 
and  a  cold  temperature  as  shown  in  Fig.  C-l.  The  remaining  box  walls  are 
perfectly  insulated.  The  objective  is  to  numerically  solve  for  the  fluid  flow 
and  heat  transfer  inside  the  box  arising  due  to  buoyancy  effects.  The 
governing  equations  for  the  steady-state  problem  assuming  laminar  flow, 
constant  properties  for  the  fluid,  no  viscous  dissipation  and  the  Boussinesq 
approximation  to  be  true  are: 


Continuity 


Du  j)v  dvv  _  n 
dx  +  3y  +  dz 


(C-l) 


Momentum  (x  direction) 


d(uu)  d(vu)  d(w  u) 
f)x  +  dy  dz 


Pr(^-  + 
dx2 


d2u  d2u  > 

dy2  dz2 


Momentum  (y  direction) 


d(U  V)  d(V  V)  d(W  v) 
dx  dy  dZ 


(C-2) 


0p  +  Pr(^+  ^4+  ^v-)+  RaPre 


dy 


dx*1  dy4 


dz* 


(C-3) 
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Adiabatic  Wall  (z=0) 


Adiabatic  Wall  (z=2) 


Fig.  C-1 .  Natural  convection  in  a  rectangular  enclosure. 


Momentum  (z  direction) 


3(uw)  0(VV\^  0{w  w) 

ax  ay  +  0z  ~ 


3P  +  Pr(^?+ 


3z 


0x^  3y‘ 


dz‘ 


(C-4) 


Energy 

3(u0)  0(_v0)  0(w0) 

Ox  +  3y  +  3z 

a2e  a2e  a2e 

ax2  ay2  az2  (c-s> 

The  above  non-dimensional  equations  have  been  derived  using  the  scaling 
of  H  and  o</H  for  the  lenght  and  velocity,  respectively;  where  H  is  the 
characteristic  dimension  of  the  box  in  the  x  direction  and  a  is  the  thermal 
difiusivity  of  the  fluid.  Pr=v/a  is  the  Prandtl  number  of  the  fluid  and 
Ra=gP(Th-Tc)H-Vav  is  the  Rayleigh  number;  v  is  the  kinematic  viscosity,  g  is 
the  gravitational  acceleration,  P  is  the  coefficient  of  thermal  expansion  and 
Th  and  Tc  are  the  temperatures  of  the  hot  and  cold  walls,  respectively.  The 
temperature  has  been  non-dimensionalized  as  0=(T-Tc)/(Th-Tc). 

As  mentioned  earlier,  the  equations  to  be  solved  are  cast  into  the 
canonical  form  in  Eq.  (1).  The  appropriate  definitions  for  p,  p,  P,  Sp  and  Sc 

are  given  in  Table  C-l: 
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Table  C-l  Definition  of  variables  in  the  canonical  equation  (Eq.  1)U) 


Equation^2) 

r 

Sp 

S  (3) 

x  momentum 

u 

Pr 

0 

0 

y  momentum 

v 

Pr 

0 

PrRa9 

z  momentum 

w 

Pr 

0 

0 

Energy 

9 

1 

0 

0 

(1)  Value  of  p  to  be  used  in  the  canonical  form  is  1,  and  is  set  in  variable 
RHOCON  in  the  BLOCK  DATA. 

(2)  It  is  not  required  to  cast  the  continuity  equation  in  the  canonical  form. 

(3)  The  pressure  gradient  terms  are  incorporated  internally  in  the 
program  code  and  need  not  be  specified  as  source  terms  in  the  subroutine 
USE. 


Taking  advantage  of  the  symmetry  plane  z=l,  the  computational  domain  is 
restricted  to  1.  The  appropriate  boundary  conditions  then  become: 


09  n 
u,v,w’  0z  =  0 


at  z=0 


0U  0v  09 
0Z’0Z’W’0Z 


=  0 


at  z=l 


u=v=w=0,  9=1  at  x=0 
u,v,w,0=O  at  x=l 
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ae  n 

u-v-w-3y  =0 


at  y=0  and  1 


Program  Implementation; 

As  mentioned  earlier,  the  user  needs  set  up  only  the  subroutine  USE  and 
BLOCK  DATA  for  a  particular  problem  to  be  solved.  The  Subroutine  USE 
and  other  BLOCK  DATA  options  to  be  used  for  this  problem  are  given  at  the 
end  of  the  source  code  in  Appendix  B.  The  subroutine  is  commented 
appropriately  for  the  user  to  understand  the  program  usage.  Special 
considerations  in  the  subroutine  USE  are  as  follows: 

(1)  Whenever  a  gradient  of  a  dependent  variable  is  specified  to  be  zero  at 
a  boundary,  the  boundary  value  of  T  associated  with  that  §  is  set  to  zero  in 
ENTRY  DIFFUS.  Thus  the  dependency  of  the  boundary  value  of  the 
variable  on  the  solution  of  the  variable  at  the  adjacent  node  is  dropped, 
leading  to  much  faster  convergence.  The  same  solution  is  obtained  if  the 
boundary  T  values  were  kept  as  those  in  Table  C-l,  however,  convergence 
will  be  slower. 

(2)  For  t  y  momentum  equation,  the  full  source  term  RaPr9  is  not 
used  right  from  the  start  of  the  iterations,  but  only  a  portion  of  it  is  used  and 
the  source  term  is  gradually  increased  with  iterations  to  its  fullest  value. 
This  helped  in  obtaining  convergence,  which  could  not  be  obtained  earlier 
with  the  full  source  term  for  all  iterations.  The  user  is  referred  to  Patankar 
[1],  Ch.  7,  for  a  detailed  discussion  on  convergence  and  underrelaxation 
techniques  for  the  source  terms. 

(3)  Since  the  u  velocity  U(I,J,K)  is  evaluated  on  the  staggered  grid  at  the 
grid-point  location  of  the  interface  of  the  control  volumes  around 
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(X(I),Y(J),Z(K))  and  (X(I),Y(J-1),Z(K)),  in  the  source  term  RaPrO,  the 
temperature  9  must  be  evaluated  at  the  interface  via  interpolation. 

Results: 

Shown  in  Fig.  C-2  is  the  local  Nusselt  number  at  the  hot  wall 
compared  with  the  solution  from  Mallinson  and  Davis  [2],  The  local 
Nusselt  number  has  been  defined  as: 

Nu(z)  =  J|(0,y,z)dy  (C.6) 


As  can  be  seen,  the  agreement  is  quite  good.  The  small  disagreement  is 
attributed  to  the  difference  in  the  numerical  procedure,  difference  in 
discretization  procedure  and  difference  in  numerical  evaluation  of  the 
Nusselt  number. 
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sse 


0.0  0.2  0.4  0.6  0.8  1.0 

z/ZL 


Fig.  C-2.  Comparison  of  local  Nusselt  numbers. 
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